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ABSTRACT 


Recent  development  of  several  unsteady  supersonic  methods 
for  computations  of  airloads  for  elastic  bodies  of  revolution, 
asymmetric  bodies  and  body-wing  configurations  are  reported. 

These  methods  include  the  Harmonic  Potential  Panel  (HPP)  method, 
the  Bundle  Triplet  Method  (BTM)  and  the  combined  method  of  BTM 
and  the  Harmonic  Gradient  Method  (HGM)  for  body-wing 
combinations.  All  methods  are  based  on  the  generic  Harmonic- 
Gradient  (H-G)  model,  which  is  essential  in  providing  accurate 
solutions  in  the  full  frequency  domain  and  the  low  Mach  number 
range . 

Extensive  comparisons  of  computed  results  obtained  from 
these  methods  show  good  correlations  with  existing  data. 
Comparison  examples  range  from  simple  cones  and  ogive  bodies  to 
Saturn  SA-1  configuration,  to  the  cylindrical  panel  membrance  and 
to  the  NACA  wing-body  combinations.  Cases  computed  yield  steady 
and  unsteady  pressures,  generalized  forces,  stability 
derivatives,  aerodynamic  dampings  and  divergence  and  flutter 
boundaries  for  these  configurations. 

The  developed  methods  have  been  validated  with  existing 
theories  or  measured  data.  For  supersonic  aeroelastic  analysis, 
these  methods  yield  results  that  are  accurate  and  cost-effective, 
thus  rendering  them  very  favorable  for  technology  transfer  and 
industry  applications. 


TABI.fi  OF  CONTENTS 


e 


LIST  OF  FIGURES  .  vi 

NOMENCLATURE  .  xi 

CHAPTER 

1.  INTRODUCTION  .  1 

1.1  Survey  of  Literature  .  2 

1 . 2  Out  line  .  3 

2.  FORMULATION  .  6 

2.1  Wind-Fixed  Coordinate  System  .  7 

2.1.1  Boundary  Conditions  .  8 

2.1.2  Pressure  Coefficient  .  9 

2.2  Body-Fixed  Coordinate  System  .  11 

2.2.1  Boundary  Conditions  .  13 

2.2.2  Pressure  Coefficient  .  15 

2.3  Pseudo-Wind-Fixed  Coordinate  System  .  16 

2.4  Generalized  Forces  and  Stability  Derivatives  .  19 

3.  STEADY  MEAN  FLOW  .  21 

3.1  Solution  of  the  Linear  Equation  .  22 

3.2  The  Nonlinear  Equation  and  Its  Solution  .  23 

3.3  Results  and  Discussions  .  26 

4.  UNSTEADY  FLOW  COMPUTATIONS  .  28 

4.1  The  Integral  Solution  .  28 

4.2  Harmonic  Gradient  Model  .  30 

4.3  Evaluation  of  Velocities  .  32 

4.4  Methods  of  Solution  .  32 

4.5  Results  and  Discussions  .  34 


4.5.1  Results  According  to  Various 

Coordinate  Systems  .  34 

4.5.2  Effects  of  Frequency  .  36 

4.5.3  Nonlinear  Results  .  38 

4.5.4  Elastic  Bodies  .  39 

5.  FLUTTER  COMPUTATION  .  41 

5.1  Flutter  Equations  .  41 

5.2  Flutter  Results  .  44 

6.  ASYMMETRIC  BODIES  .  46 

6.1  Formulation  .  46 

6.1.1  Boundary  Conditions  .  47 

6.1.2  Pressure  Coefficient  .  49 

6.2  Mean  Flow  Solution  .  51 

6.3  Unsteady  Flow  Solution  .  58 

7.  Bundled  Triplet  Method  .  64 

7.1  New  Development  of  BTM  .  64 

7.2  Formulation  .  64 

7.3  Harmonic  Gradient  Method  .  66 

7.4  Least  Square  Procedure  .  69 

7.5  Panel  Flutter  .  72 

7.6  Results  and  Discussion  .  74 

7.6.1  Asymmetric  Bodies  .  74 

7.6.2  Cylindrical  Panel  Flutter  .  76 

7.6.3  Salient  Feature  of  BTM  .  76 

8.  Body-wing  Combinations  .  78 

8.1  Unsteady  Lifting  Surface  Method  .  78 

8.2  Body-Wing  Formulation  .  79 

i  v 


8.3  Pressure  Coefficients  .  80 

8.3.1  Unsteady  Pressure  Coefficient  for 

the  Body  .  80 

8.3.2  Unsteady  Pressure  Coefficient  for 

the  Wing  .  81 

8.4  Generalized  Forces  and  Stability  Derivatives  .  82 

8.5  Results  and  Discussion  .  83 

9.  Conclusion  . 84 

REFERENCES  .  120 

APPENDIX  .  88 

A.  LINEARIZED  EQUATIONS  IN  THE  BODY-FIXED 

COORDINATE  SYSTEM  .  88 

B.  EVALUATION  OF  THE  UNSTEADY  POTENTIAL  AND 

VELOCITIES  APPLYING  THE  HARMONIC  GRADIENT  MODEL  .  95 

C.  ON  THE  APEX  SINGULARITY  FOR  OSCILLATING  POINTED 

BODIES:  WIND-FIXED  VERSUS  BODY-FIXED  COORDINATE 

SYSTEMS  .  104 

D.  EVALUATION  OF  THE  INDUCED  UNSTEADY  POTENTIAL 

AND  VELOCITIES  BY  A  DISTRIBUTION  OF  SOURCES  .  114 

E.  EXPLICIT  EXPRESSION  OF  MATRICES  [A];  [y]^ 

AND  [H]  t  IN  EQ.  (7.11)  .  118 


v 


LIST  OF  FIGURES 


Figure  Page 

1.  Body  bending  flutter  of  British  J.T.V.I.  ramjet 

missile  .  125 

2.  Sequence  of  flutter  modes  at  intervals  of  1/100 

second  during  the  flutter  incident  .  126 

3.  Wind-fixed,  body-fixed,  pseudo-wind-fixed  coordinate 

system  .  127 

4.  Comparison  of  normalized  velocities  on  a  cone-cylinder 

surface  at  =  2.075  and  angle  of  attack  a=0°  .  128 

5.  Mean  flow  pressures  for  a  parabolic-ogive  at  M^  =  2.0 

and  angle  of  attack  o  =  0°  129 

6.  Mean  flow  pressures  for  a  parabolic-ogive  at  M^  =  3.0 

and  angle  of  attack  a  =  0°  .  130 

7.  Total  pressure  distributions  of  an  ogive-cylinder- 

boattail  body  at  M^  =  3.0  .  131 

8.  Panel  arrangement  for  axisymmetric  bodies  .  132 

9.  Harmonic  gradient  model  .  133 

10.  In-phase  and  out-of-phase  pressure  coefficients  for 
a  cone  in  pitching  mode  at  =  2.0  and  reduced 

frequency  k  =  2.0  .  134 

11.  In-phase  and  out-of-phase  pressure  coefficients  for 

a  cone  in  first  bending  mode  at  Mw  =  2.0  and  a  reduced 

frequency  k  =  2.0  .  135 

12.  In-phase  and  out-of-phase  pressure  coefficients  for  a 
cone  in  second  bending  mode  at  M^  =  2.0  and  reduced 

frequency  k  =  2.0  .  136 

13.  Modulus  and  argument  of  the  generalized  forces  versus 

Mach  number  for  a  cone  in  first  bending  mode  (1=3)  at 
reduced  frequency  k  =  2.0  .  137 

14.  Modulus  and  argument  of  the  generalized  forces  versus 

Mach  number  for  a  cone  in  first  bending  mode  (1=4)  at 
reduced  frequency  k  =  2.0  .  138 


Figure 


Page 


15.  Modulus  and  argument  of  the  generalized  forces  versus 
Mach  number  for  a  pa r ab o 1  i c-og i ve  in  first  bending  mode 

(1  =  3)  at  reduced  frequency  k  =  2.0  .  139 

16.  Modulus  and  argument  of  the  generalized  forces  versus 
Mach  number  for  a  parabol ic=ogive  in  second  bending 

mode  (1=4)  at  reduced  frequency  k  =  2.0  .  140 

17.  Modulus  and  argument  of  the  generalized  forces  versus 
Mach  number  for  a  cone-cylinder  in  first  bending  mode 

(1  =  3)  at  reduced  frequency  k  =  2.0  .  141 

18.  Modulus  and  argument  of  the  generalized  forces  versus 
Mach  number  for  a  cone-cylinder  in  second  bending 

mode  (1=4)  at  reduced  frequency  k  =  2.0  .  142 

19.  Comparison  of  theoretical  damping-in-pitch  moment 

coefficients  for  a  parabolic-ogive  at  various  Mach 
numbers  .  143 

20.  Comparison  of  theoretical  and  experimental  damping- 

in-pitch  moment  coefficients  for  a  parabolic-ogive 
cylinder  at  various  Mach  numbers  .  144 

21.  Comparison  of  theoretical  and  experimental  damping- 

in-pitch  moment  coefficients  for  a  cone  frustum  at 
various  Mach  numbers  .  145 

22.  Modulus  and  argument  of  the  generalized  forces  versus 
reduced  frequency  for  a  cone  and  a  parabolic-ogive 

in  pitching  mode  at  M^  =  2.5  .  146 

23.  In-phase  and  out-of-phase  pressure  coefficients  for  a 
parabolic  ogive  in  pitching  mode  at  =  1.5  and  high 

reduced  frequencies  .  147 

24.  Comparison  of  the  in-phase  and  out-of-phase  pressure 
coefficients  at  M  =  1.5  for  a  5.7°  cone  and  a  flat 

90 

plate  pitching  at  the  apex  . 148 

25.  D amp i n g- i n -p i t ch  normal  force  coefficients  versus 

Mach  number  for  a  cone  and  a  parabolic-ogive  .  149 

26.  D amp i n g- i n-p i t ch  moment  coefficient  for  a  parabolic- 

ogive-cylinder  at  various  Mach  numbers  .  150 

27.  Aerodynamic  damping  coefficient  versus  Mach  number 

for  a  cone  cylinder  .  151 


Figure 


28.  Pressure  coefficients  for  a  Saturn  SA-1  configuration 
at  =  2.0,  reduced  frequency  k  =  1.8  and  center  of 

rotation  at  the  apex  .  152 

29.  Comparison  of  computed  and  experimental  aerodynamic 
damping  coefficients  versus  Mach  number  for  a  Saturn 

SA-1  configuration  vibrating  in  first  bending  mode  ....  152 

30.  Schematic  view  of  body  in  pitching  and  plunging  motion  154 

31.  Flutter  boundaries  for  a  7.5°  cone  at  M  =  2.0  .  155 


32.  Flutter  boundaries  for  a  7.5°  cone  at  M 

oo 


3 . 0 


156 


33.  Flutter  speed  boundaries  versus  Mach  number  for 

a  7.5°  cone  . 

34.  Flutter  frequency  boundaries  versus  Mach  number 

for  a  7.5°  cone  . 

35.  Divergence  boundaries  versus  Mach  number  for  a 

7.5°  cone  . 


157 


158 


159 


36.  Panel  and  control  points  arrangement  for  asymmetric 

bodies  .  160 

37.  Pressure  coefficient  for  an  asymmetric  cone  at  M^  =  2.0, 

a  =  0°  .  161 

38.  Pressure  coefficient  for  a  half  circle;  half  ellipse 

cone  at  M^  =  2.0,  a  =  0°  .  162 


39.  Normal-force  coefficient  for  elliptic  cones  versus  the 

ellipticity  ratio  a/b  at  M^  =  3.0  .  163 


40.  Pitching  moment  coefficient  for  elliptic  cones  versus 

the  ellipticity  ratio  a/b  at  M^  =  3.0  .  164 

41.  D amp i n g- i n-p i t ch  normal  force  coefficient  and  damping- 
in-pitch  moment  coefficient  for  elliptic  cones  versus 

the  ellipticity  ratio  a/b  at  =  3.0  .  165 

42.  In-phase  and  out-of-phase  pressure  coefficient  for 
elliptic  cones  in  pitching  mode  =  3.0  at  and  reduced 

frequency  k  =  1.0  along  the  x-axis  .  166 


v  i  i  i 


Figure 


Page 


43.  In-phase  and  out-of-phase  pressure  coefficient  for 
elliptic  cones  in  first  bending  mode  at  =  3.0  and 

reduced  frequency  =  1.0  along  the  x-axis  .  167 

44.  In-phase  and  out-of-phase  pressure  coefficient  for 
elliptic  cones  in  second  bending  mode  at  =  3.0  and 

reduced  frequency  =  l.u  along  the  x-axis  .  168 

45.  In-phase  and  out-of-phase  pressure  coefficient  for 
elliptic  cones  in  pitching  at  mode  =  3.0  and  reduced 

frequency  k  =  1.0  along  the  polar  angle  8  .  169 

46.  In-phase  and  out-of-phase  pressure  coefficient  for 
elliptic  cones  in  first  bending  mode  at  =  3.0  and 

reduced  frequency  k  =  1.0  along  the  polar  angle  8  .  170 

47.  In-phase  and  out-of-phase  pressure  coefficient  for 
elliptic  cones  in  second  bending  mode  at  Mm  =  3.0  and 

reduced  frequency  k  =  1.0  along  the  polar  angle  8  .  171 

48.  Pressure  distribution  per  unit  angle  of  attack  along 
the  y-axis  for  a  circular  cone,  elliptic  cone  and  a 

flat  plate  at  =  2.0  .  172 

49.  Sketch  showing  bundled  triplet  arrangements  .  173 

50.  Sketch  showing  the  domain  of  influence  based  on  the 

bundled  triplet  method  .  174 

51.  Steady  pressure  distributions  for  an  asymmetric  cone 

at  =  2.0  and  angle  of  attack  a  =  0°  175 

52.  Steady  pressure  distributions  for  an  asymmetric  cone 

at  =  2.0  and  angle  of  attack  a  =  0°  176 

53.  Steady  pressure  distributions  for  an  elliptic  cone 

at  =  2.0  and  angle  of  attack  a  =  0°  177 

54.  Stability  derivatives  of  a  family  of  elliptic  cones 

at  various  Mach  numbers  .  178 

55.  Cylindrical  flutter  showing  two  vibrating  modes  .  179 

56.  Real  part  of  the  generalized  aerodynamic  force  Qjj 

versus  circumferential  mode  number  .  180 


Figure  Page 

57.  Imaginary  part  of  the  generalized  aerodynamic  force 

Qj j  versus  circumferential  mode  number  .  181 

58.  Sketches  of  wing-body  configurations 

a)  Aspect  ratio  AR  =  3  triangular  wing-body 

b)  Aspect  ratio  AR  =  3  swept  wing-body  .  182 

59.  Static  and  dynamic  moment  derivatives  for  an  aspect 

ratio  AR  =  3  triangular  wing-body  .  183 

60.  Static  and  dynamic  moment  derivatives  for  an  aspect 

ratio  AR  =  3  swept  wing-body  .  184 

61.  Conical  coordinate  and  cylindrical  coordinate  systems 

for  a  circular  cone  . .  185 


x 


NOMENCLATURE 


S  yrobo  1 
a 


CM 


a 


c 


p 


Definition 

speed  of  sound 

pitching  moment  coefficient 

damp i n g- in-p i t ch  moment  coefficient 

P  -  P 

M 

pressure  coefficient,  \  /2p — LP" 


C° 

p 

C1 

p 

F(«) 
f  (  <) 

M«> 

g(x) 
h  ( t ) 
h  (  x  ,  t  ’  ) 
I 

a 

K 

a 

k 

L 

M 

M 


mean  flow  pressure  coefficient 

unsteady  flow  pressure  coefficient 

dipole  strength  distribution 
source  strength  distribution 
vortex  multiplets  strength  distribution 
normalized  natural  mode  shape 

generalized  coordinate  for  plunging  oscillation 
instantaneous  normal  displacement  of  the  body 
mass  moment  of  inertia  of  the  body 

stiffness  of  the  body  in  plunging 
stiffness  of  the  body  in  pitching 

w  L 

reduced  frequency,  g— 

body  length  (dimensional) 
free-stream  Mach  number 


m 


total  mass  of  the  body 


n  outward  normal  to  the  body  surface 

p ( x )  normalized  natural  mode  shape  expressed  in  the 

pseudo-wind-fixed  coordinate  system 


®I  J 

R 

r 

a 

S  =  0 

s 

a 

t 

t  ' 

u 


generalized  forces 
body  radius 

dimensionless  radius  of  gyration  about  pitch 

axis,  1  *  a 
mL2 

body  surface 

static  moment  of  the  body 


time  measured  in  a  body-fixed  reference  frame 
time  measured  in  a  wind-fixed  reference  frame 
f ree-s t ream  speed  (dimensional) 


V 

-¥ 

VB 

XG 


x,  r,  B 


x,  y,  z 
x'  ,  r’  ,  8  ‘ 
x’  ,y'  ,z‘ 


total  velocity  vector 

velocity  vector  of  the  body  surface 
pitching  axis  location  of  the  body 
body-fixed  cylindrical  coordinates 
body-fixed  cartesian  coordinates 
wind-fixed  cylindrical  coordinates 
wind-fixed  cartesian  coordinates 


Greek  Symbols 

<*(t)  generalized  coordinate  for  pitching  oscillation 

p  /M2  -  T 

oo 


f  ratio  of  specific  heats 

5 ( t ’ )  instantaneous  amplitude  of  oscillation 

S 0  amplitude  of  oscillation 

x  flexible  part  of  the  unsteady  flow  potential  in 

the  pseudo-wind-fixed  coordinate  system 

f  source  of  dipole  coordinate 

free-stream  density 


X  1  1 


$  complete  time-dependent  perturbation  potential 

mean  flow  perturbation  potential 
mena  flow  first-order  iteration  potential 
cross  flow  perturbation  potential 

X0  homogeneous  solution  for  the  mean  flow  second- 

order  iteration  equation 

*  rigid  part  of  the  unsteady  flow  potential  in  the 

pseudo-wind-fixed  coordinate  system 

f0  particular  solution  to  the  mean  flow  second-order 

iteration  equation 

Q  full  potential 

w  frequency  of  oscillation  (dimensional) 

wh  uncoupled  angular  frequency  in  vertical 

translation 

uncoupled  angular  pitching  frequency 

V  gradient  operator 

?2  Laplacian  operator 

V*  divergence  operator 

Subscripts  and  Superscripts 

()',()", etc .  derivatives  with  respect  to  the  independent 
variable 

( ) x , r , g , t  partial  derivatives 

(  •) 


defivative  with  respect  to  time 


CHAPTER  1 


INTRODUCTION 


Long  slender  missiles  or  rockets  at  cruising  supersonic 
speeds  are  susceptible  to  a  number  of  aeroelastic  instabilities. 
First,  it  is  known  that  the  stability  and  control  characteristics 
of  high  speed  flexible  bodies  may  be  significantly  influenced  by 
the  distortion  of  the  structure  under  transient  loading 
conditions.  Second,  the  body/fin  configurations  are  likely  to 
flutter  as  a  result  of  the  properly  phased  short-period  rigid 
body-fin  mode  and  the  body  bending  mode.  In  fact,  the  earlier 
flutter  incident  of  the  British  JTVI  ramjet  missile  in  1955,  as 
shown  in  Figs.  1  and  2,  is  one  such  example  (Ref.  1). 

Furthermore,  the  problem  of  the  store-airframe  interaction, 
during  the  cruise  and/or  maneuver  phase,  of  modern  aircraft  has 
been  a  major  concern  for  design  and  performance.  The  effects  of 
this  type  of  interaction  could  sometimes  change  the  airload,  and 
hence  the  wing  flutter  characteristics,  rather  drastically.  For 
example,  problems  such  as  stores  in  pitch-yaw  combined 
oscillations  and  the  tip-missile  influence  are  among  the  critical 
factors  related  to  aircraft  flutter.  Clearly,  the  prediction  of 
these  boundaries  relies  almost  exclusively  on  the  unsteady 
aerodynamic  inputs.  The  objective  of  the  present  work  is 
therefore  to  provide  a  generalized  Harmonic  Potential  Panel  (HPP) 
method  for  computing  the  unsteady  aerodynamics  of  arbitrary 
flexible  bodies  and  body-fin  configurations  in  the  supersonic 
flow  regime. 
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1.1  Survey  of  Literature 


Currently,  several  panel  methods  have  claimed  success  for 
predictions  of  the  steady  aerodynamics  for  wing-body  combinations 
in  the  supersonic  flow  regime  (Refs.  2  and  3).  In  recent  years, 
computational  methods  for  unsteady  supersonic  flow  prediction 
have  been  extensively  investigated  (Ref.  4).  However,  these 
approaches  are  mostly  formulated  for  wing  planform  calculations. 
On  the  other  hand,  the  investigation  of  unsteady  supersonic  flow 
for  oscillating  bodies  in  the  past  has  been  mostly  based  on 
slender  body  or  not-so-s lender  body  theories  for  rigid-body 
oscillations  in  the  low-frequency  range. 

In  the  slender  body  limit  the  Adams  and  Sears  theory  (Ref. 

5)  was  extended  to  unsteady  flow  by  Garrick  (Ref.  6).  Although 
Garrick’s  theory  is  valid  for  all  frequencies  and  for  flexible 
bodies,  it  has  the  limitations  of  being  independent  of  the  Mach 
number  and  it  is  found  too  inaccurate  for  bodies  of  practical 
thickness . 

In  the  not-so-s 1 ender  body  limit,  Lansing  (Ref.  7)  used  a 
frequency  expansion  procedure  for  treatment  of  rigid-body 
oscillations,  and  Platzer  and  Sherer  (Ref.  8)  applied  the 
linearized  method  of  characteristics  (LMOC)  for  rigid  bodies  in 
low-frequency  oscillations.  Tobak  and  Wehrend  (Ref.  9)  extended 
Van  Dyke’s  first-  and  second-order  theories  (Ref.  10)  to  unsteady 
flow  for  cones.  Bond  and  Packard’s  theory  (Ref.  11)  for  flexible 
bodies  appeared  in  1961;  however,  it  was  found  to  involve 
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erroneous  boundary  conditions,  as  was  pointed  out  by  Hoffman  and 
Platzer  ( Ref.  12 ) . 

It  appears  that  little  progress  has  been  made  in  the 
development  of  a  viable  computational  method  for  wing-body  or 
body-fin  combinations,  for  unsteady  supersonic  aerodynamic 
predictions,  which  is  uniformly  valid  in  the  complete  frequency 
domain.  Hence,  the  present  work  consists  of  the  development  of  a 
viable  method  in  the  f u 1 1- f r equency  domain  for  computations  of 
the  unsteady  aerodynamics  of  arbitrary  flexible  bodies  in  the 
supersonic  flow  regime. 


1.2  Outline 

For  bodies  in  pitch  motion,  the  proper  choice  of  the 
coordinate  system  has  been  subject  to  some  controversy  in  the 
past.  In  Chapter  2  the  formulations  in  the  wind-fixed,  the  body- 
fixed  (Appendix  A),  and  the  pseudo-wind-fixed  coordinate  systems 
are  presented. 

In  Chapter  3  the  solution  procedure  for  the  steady  mean  flow 
problem  is  developed.  The  linearized  equation  for  the  mean  flow 
is  solved  by  using  the  Karman  and  Moore  procedure  (Ref.  13).  The 
nonlinear  equation  for  the  mean  flow  is  solved  according  to  Van 
Dyke’s  iterative  scheme  (Ref.  10).  Computed  mean— flow  velocities 
and  pressures  for  a  cone-cylinder,  a  parabolic-ogive,  and  a 
par abo  1  i c-og i ve-boa 1 1 a i 1  bodies  are  presented. 

In  Chapter  4  the  solutions  of  the  equations  in  various 
coordinate  systems  as  derived  in  Chapter  2  are  obtained.  The 
Harmonic  Gradient  concept  is  applied  to  the  dipole  strength  so 
that  the  number  of  panel  elements  becomes  least  affected  by  the 
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given  Mach  number  and  reduced  frequency.  It  is  shown  that  when 
the  wind-fixed  system  is  applied  in  a  straight-forward  manner, 
the  solution  obtained  is  totally  contaminated  by  a  spurious 
leading-edge  singularity.  In  Appendix  C,  it  is  shown  that  such 
singularity  is  totally  removed  in  a  conical  coordinate  system. 
Effects  of  frequency  and  mean-flow  nonlinearities  are 
investigated  in  the  body-fixed  coordinate  system.  Comparisons  of 
the  present  results  with  NASA’S  measured  data  (Refs.  14  and  15) 
of  the  aerodynamic  damping  coefficients  for  different  bodies  are 
shown . 

In  Chapter  5  the  equations  for  flutter  computations  are 
presented.  Flutter  boundaries  for  a  7.5*  cone  are  determined  by 
using  the  linear  and  nonlinear  method.  Comparisons  of  the 
present  results  with  NASA’s  measured  data  (Ref.  16)  are 
presented . 

In  Chapter  6  the  formulation  for  bodies  with  asymmetric 
cross-section  is  presented  in  the  body-fixed  coordinate  system. 
For  ease  of  application  a  spline  panel  method  is  used  to 
determine  the  pressure  and  forces  acting  on  the  body. 

Comparisons  of  the  mean  flow  pressures  and  static  forces  for 
elliptic  cones  with  those  computed  by  USSAERO  code  (Ref.  2)  show 
clear  limitations  of  the  present  method  when  the  asymmetry  in  the 
body  cross-section  increases. 

In  Chapter  7,  the  development  of  the  Bundle  Triplet  Method 
(BTM)  is  presented.  The  BTM  is  a  more  general  method  in 
treatments  of  asymmetric  bodies  in  that  the  line  doublet 
formulation  of  the  spline  panel  method  is  generated  to  a 
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generalized  triplet  one.  The  body  cross-section  is  also  divided 
into  the  so-called  ”pie  sectors"  but  a  line  source  and  a  line 
doublet  are  superposed  in  each  sectors;  thus  forming  a  bundle  of 
triplet  lines.  A  least  square  procedure  is  used  to  account  for 
the  interference  between  sectors.  Substantial  improvements  are 
found  in  the  results  over  those  obtained  using  the  methods  in 
Chapter  6.  Moreover,  the  BTM  is  validated  through  the 
application  cases  of  cylindrical  panel  flutter  in  which  the 
generalized  forces  obtained  are  in  excellent  agreement  with  exact 
theories  in  the  higher  frequency  range  and  for  higher  order 
modes . 

In  Chapter  8,  the  formulation  for  the  body-wing  combinations 
is  presented.  The  acceleration  potential  version  of  the  Harmonic 
Gradient  Method  (AHGM)  is  in  combined  use  for  computations  of  the 
body-wing  aerodynamics.  The  present  formulation  not  only  allows 
the  computed  domain  to  be  confined  to  the  surface  panels  but  it 
also  directly  yields  the  unsteady  pressure  on  the  wing  surface. 
Therefore,  the  present  method  is  a  very  effective  one  in 
treatments  of  unsteady  aerodynamics  over  body-wing 
configurations. 

Finally,  conclusions  from  the  present  investigation  and 
recommendations  for  future  work  are  presented  in  Chapter  9. 

In  what  follows,  all  variables  are  nond imens i ona 1 i z ed  by  the 
true  length,  and  time  scales  defined  by  the  body  length,  and  the 
body  length  divided  by  the  freestream  velocity,  respectively. 
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CHAPTER  2 


FORMULATION 

In  this  chapter,  the  formulation  for  a  body  of  revolution 
performing  bending  oscillations  is  presented. 

The  fluid  flow  is  assumed  to  be  inviscid  and  isentropic. 
Thus,  the  fluid  velocity  V  can  be  defined  by  the  scalar  potential 
Q( x ’ , y ’  , z ’  , t ’  )  ,  i.e. 

V  =  70  (2.1) 

The  governing  equation  for  the  potential  0  is  the  full 
potential  equation 

a2  72  Q  =  Qt,t,  +  |p-  (70)2  +  ( 70  •  7)  ( 70 ) 2  (2.2) 

where  "a"  is  the  local  speed  of  sound  given  by 

a2  =  jL-  -  (f-1)  {ot,  +  \  [ ( 70) 2  -  1]}  (2.3) 

oo 

The  time  coordinate  above,  t’,  is  based  on  a  spatial 
coordinate  system,  (x’,y’,z’),  that  is  fixed  with  respect  to  the 
fluid  at  infinity. 

It  is  known  that  various  linearized  small  disturbance 
equations  can  be  derived  from  Eqs.  (2.2)  and  (2.3)  depending  on 
the  coordinate  system  chosen.  Their  formulations  according  to 
the  wind-fixed,  the  body-fixed  and  the  pseudo-wind-fixed 
coordinate  systems  are  discussed  in  order. 
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2.1  Wind-Fixed  Coordinate  System 
If  a  cylindrical  coordinate  system  with  its  x’-axis  parallel 
to  the  freestream  velocity  is  adopted  (see  Fig.  3)  then  the  full 
potential  can  be  written  as 

Q  (x’ , r’ ,©’ , t’ )  =  x’  +  4>  (x’.r’.O’.t’)  (2.4) 


where  <S  (  x  ’  ,  r  ’  ,  © 1  ,  t  ’  )  is  the  perturbation  potential. 

Substituting  £q.  (2.4)  into  Eqs.  (2.2)  and  (2.3)  and 
retaining  only  linear  terras  in  4  yields  the  following  linear 
equat i on 


(1  -  Mf) 


+  *r 
r  r 


-  2M2 


-  =  0  (2.5) 


For  a  body  of  revolution  in  unsteady  motion  with  small 
amplitude  of  oscillation  the  perturbation  potential  can  be  split 
into  two  parts;  the  mean  flow,  zero-angle-of-attack ,  potential, 
0o(x’,r’),  and  the  unsteady  flow  potential,  # t ( x ’ , r ’ ) cos© ’ 

♦  (x  ’  ,  r  ’  ,  ©’  ,  t  ’  )  =  0O  (  x  ’  ,  r  ’  )  +  <5  ( t  ’  )*x  (x’  ,  r’  )  cos©’  (2.6) 


where  <5  ( t  ’  )  is  the  instantaneous  amplitude  of  oscillation.  For 
simple  harmonic  oscillations  5 ( t  ’  ) = S0 e 1 k 1  ’  ,  where  k  is  the 
reduced  frequency.  Substituting  these  relations  into  Eq.  (2.5) 
and  collecting  terras  of  order  one  and  50 ,  one  obtains 


t 

\ 


1 


M*)  *o: 


+  7>*0r 


0 


(2.7) 
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(  1  '  ML)  ^1  X  *  X  '  +  ^1  r  •  r  •  +  fT  ^lr- 

~  \tt  * i  -  2ikM2  *lx,  +  M^k2  ^ i  =  0  (2.8) 


These  are  the  governing  equations  in  the  wind-fixed  coordinate 
system  for  <fiQ  and  <j>  x  . 


2.1.1  Boundary  Conditions 

Physical  considerations  suggest  that  all  perturbations 
should  vanish  upstream  of  the  Mach  wave  emanating  from  the  body 
apex  , 


*0  (  x  ’  <  r  ’  ) 
*x  (x*  ,  r’  ) 


•P0%,  (x’.f’)  =  *0  r  ,  (  x  ’  ,  r  *  ) 
^lx.  (x  ’  ,  r  ’  )  =  <j>x  r  ,  (x  ’  ,  r  ’  ) 


0 


at  x’-/Jr’<0 


0 


(2.9) 


and  that  at  any  instant  the  flow  must  be  tangent  to  the  body 
surface.  In  the  wind-fixed  coordinate  system  the  latter 
condition  can  be  expressed  as 

=0  at  S  =  0  (2.10) 


where  S ( x ’ , r ’ , © ’ , t ’ ) =0  describes  the  body  surface.  For  a  body  of 

revolution  performing  small  amplitude  bending  oscillations,  the 

equation  of  the  body  surface  can  be  simplified  (Ref.  17)  to 

S(x’  ,r’  ,9’  ,  t’)  =  r’-R(x’)  +  5 ( t ’ ) g ( x ’ ) cos  9 ’  +  0 ( Sg  )  =  0 

(2.11) 

where  g(x’)  is  the  normalized  natural  mode  shape. 
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Now,  substituting  Eqs.  (2.11),  (2.4)  and  (2.6)  into  (2.10) 

yields 

*0r •  '  R’  (1  +  *ox  ■  )  +  ) cose’ 

(g’(l+*0x.)  -  R’*lx.  +  *lr.  +  ikg]  +  0(i§)  =  o  (2.12) 

at  r ’  =  R( x ’ ) -S ( t ’ ) g(x ’ ) cose’ 

It  is  now  necessary  to  remove  the  implicit  dependence  on 
5 ( t ’ )  in  the  potential  functions  by  expressing  the  boundary 
condition  at  the  mean  position  (Ref.  17)  r’=R(x’).  This  is 
accomplished  by  using  a  Tay ior-ser ies  expansion  about  r’=R(x’). 
Performing  this  expansion  in  Eq.  (2.12)  and  collecting  terms  of 
order  one  and  SQ  ,  yields 

*0r ■  ~  R’*ox -  =  H’  at  r’  =  R(x’ )  (2.13) 

r  *  "  R’^lx'  =  -g 1 (x ’ ) ( 1  +  ^Q x ,  )  +  g(x’). 

(0Or ' r '“R’^ox ’ r ’ )  “  ) 

at  r’  =  R ( x ’ )  (2.14) 

2.1.2  Pressure  Coefficient 

Based  on  the  time-dependent  Bernoulli  equation,  the  exact 
isentropic  pressure  coefficient  is  expanded  to  yield  the  mean 
flow  and  unsteady  flow  pressure  coefficients, 

Cp  =  C°  +  Cl  5 0  e1 kl ’ cos9’  (2.15) 

and  after  Tay lor-ser ies  expansion  about  the  mean  position 
r’=R(x’),  C£  and  can  be  expressed  as 
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c°p  =  —  [s/-l] 

i  M2 


at  r ’ - R ’ ( x )  (2.16) 


•2S0  {(1  +  0O*.  )  0lx>  +  ^or'^lr'  +  ik*i  _  g(x') 

If),  <  r  ■  (1+*ox  ’  )  +  *0r  •  r  '  *  0  r  '  H  at  >"  =  R(x)  (2.17) 


and 


■[ 


1  "  i—L  2^0X  •  +  *0x  ■  +  ^Or') 


1/r-l 


(2.18) 


As  shown  in  Appendix  C,  in  the  slender  body  approximation, 
>0r  is  given  by 


* 


RR’ 


0  r 


and  its  derivatives  with  respect  to  x’  and  r’  are  respectively 


R’2  +  RR’ ’ 


Or  1  « 


,  .  _  _  RR’ 

,  and  0O  r  .  r  ■  _  >  2 


At  the  body  apex  both  terms  behave  like 

R’  (0) 

RIO) 

thus,  if  R’(0)  does  not  vanish,  the  second  order  derivatives 
#0r,x,  and  <fi0r,T,  are  singular  at  the  body  apex  since  R(0)=0. 
Hence,  the  second  order  derivative  terms  in  Eqs.  (2.14)  and 
■2.17)  associated  with  g(x’)  will  result  in  an  apex  singularity 
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if  g(0)  does  not  vanish.  The  cause  of  this  apex  singularity  was 
discussed  by  Platzer  and  Liu  (Ref.  17),  where  they  pointed  out 
that  it  was  due  to  the  assumption  that  gz  ( x >  ) / R2 ( x >  ) « \  was 
not  uniformly  valid  at  the  apex  for  a  nonvanishing  g(0). 
Apparently,  this  singularity  is  inherent  to  the  thickness  part  of 
the  solution,  since  in  the  slender  body  limit,  Hoffman  and 
Platzer  (Ref.  12)  have  shown  that  the  unsteady  pressure  remains 
regular.  It  is  shown  in  Appendix  C  that  this  finding  also 
applies  to  not-so-s lender  bodies.  Hence,  to  circumvent  this 
singularity  in  the  formulation,  the  obvious  choice  is  the  body- 
fixed  coordinate  system. 

2.2  Body-Fixed  Coordinate  System 

As  shown  in  Fig.  3  the  body-fixed  coordinate  system  requires 
that  the  x-axis  remain,  at  all  times,  the  axis  of  the  body, 
whereby  each  right  cross  section  is  circular  and  contains  the  r- 
ax  i  s  . 

Let  h(x,t’)  be  the  instantaneous  normal  displacement  of  the 
body,  i  . e . 


h  (  x  ,  t  ’  )  =  <5  0  e 1  *  1  g  (  x  ) 


(2.19) 


The  full  potential  in  a  cartesian  coordinate  system  (x,y,z) 
(see  Appendix  A)  can  be  written  as: 

Q(  x , y , z , t )  =  x  +  hx(x,t)z  +  $ ( x , y , z , t )  (2.20) 


The  linearized  equation  for  $  as  derived  in  Appendix  A  is 


$ 

X  X 

2^X  X 

z*xx 

-  hx 

Z  $  4-  $ 

xx  x  yy 

4- 

hxx 

*z 

-  MiO 

t  t 

kx  t 

t  z*x 

-  hx  t  z$ x t  + 

K 

t*. 

4- 

ht  *,t 

4- 

2* 

X  t  ~ 

2hx 

x  t.  2  ^  X 

-  2hxxz*xt 

4- 

2hx 

t*2 

+ 

2h 

<t> 

X  z  t 

-  2hxt2$ 

x  x  +  2ht  4>x  2 

4- 

$ 

XX 

(1 

-  2hxxz) 

- 

hx 

X  x  Z* 

+ 

X 

h  $ 

XX  2 

+  2*xlhx] 

(2.21) 

Clearly, 

i  t 

can 

b  e 

seen 

that  the  linearized 

equation 

in  the 

fixed  coordinate  system  is  different  from  the  one  obtained  in  the 
wind-fixed  system. 

As  derived  in  Appendix  A,  Eqs.  (A. 18)  and  (A. 19),  the 
governing  equations  for  $Q  and  $x  for  a  body  of  revolution  in  the 
cylindrical  coordinates  (x,r,9)  are 


(1  -  M^)  $ 0  x  x  +  <p 0  r  r  +  ~  r  r  =  G 


(2.22) 


(1  —  M2  )  0,  4-  d> .  4-  —  <h 

v  o*7  flxx  ^lrr  r  "lr 


fT  ~  2ikMl^i x  +  Mtk20i  -  Gi(g>0o) 


(2.23) 


where  Gx  represents  the  mean  flow  and  the  flexible  mode 
interaction  and  is  given  by 


=  M2 


2^(^0r  -  («’  +  ik«) 


(^or  ~  r^o«^H«k2  +  g”  ) 

g  "  (  * 0  x  x  r  -  0O  r  ) 


+  r— (  g  ’  ’  ) 

ax'0*8 


(2.24) 
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Eq.  (2.23)  contains  Van  Dyke’s  (Ref.  10)  steady  angle  of 
attack  equation  and  Hoffman  and  Platzer’s  (Ref.  18)  low-frequency 
equation  as  special  cases.  It  should  be  pointed  out  that  Eq. 
(2.23)  for  steady  angle  of  attack  differs  from  McCanless '  (Ref. 
19)  equation  which  is  clearly  in  error  since  his  formulation 
starts  from  a  linearized  equation  rather  than  the  full  potential 
equat  ion . 

2.2.1  Boundary  Conditions 

The  boundary  conditions  at  the  apex  are  the  same  as  those 
given  by  Eq.  (2.9),  while  the  tangential  boundary  condition  in 
the  body-fixed  system  states  that  the  relative  normal  velocity  of 
a  fluid  particle  to  the  body  surface  is  zero  at  any  instant,  i.e. 

(VQ  -  VB ) *n  =  0  at  S  =  0  (2.25) 

where  VB  is  the  velocity  of  the  body  surface,  S  is  the  body 
surface  which  for  a  body  of  revolution  in  a  cylindrical 
coordinate  system  can  be  expressed  as 

S  =  r-R(x)  =  0  (2.26) 

and  n  is  the  outward  normal  to  the  body  surface.  For  a  body 
performing  bending  oscillations,  VB  can  be  shown  to  be 

VB  =  5 ( t ) [ g ’  Rcos©  ex  -  g ( cos©  er  -  sin©  e  )]  (2.27) 

The  total  velocity  in  this  curvilinear  coordinate  system  is 
obtained  (see  Appendix  A)  as 
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7Q  =  [  1  +•  0Qx  -  $(t)  g’  ’  0OxRcose  +  i(t)^lxcos9]ex 
+  [<P0r  +  8  (  t)^lr  cos©  +  g’  S  (  t  )  cos6]er 

-  <5(t)[g’  sine  +  g-  sine]eQ  (2.28) 

Now  substituting  Eqs.  (2.28)  and  (2.27)  into  (2.25)  yields 
[1  -  <5  ( t )  g  ’  ’  R  cose  +  ^ax  +•  5(t)^lxcos6 
-  i(t)g’  R  cose]  ( -  R  ’  )  +  *>0r 

+  (5(.t)g*  +  5(t)^lr  +  <5  ( t )  g )  cos  e  =  0  at  r  =  R(x) 

(2.29) 


where  for  simple  harmonic  oscillations  4(t)=ikS(t)  and 
<$  =  J0eikt.  Substituting  this  relation  into  Eq.  (2.29)  and 
collecting  like  terms  of  order  one  and  <5 0  ,  yield 

*or  -  R’*ox  =  R’  (2.30) 

*lr  -  R’*ix  =  -  +  g’RR’)  -  g”HR’^0x 

at  r=R ( x )  (2.31) 


Eq.  (2.31)  reduces,  to  the  Lighthill’s  (Ref.  20)  boundary 
conditions  when  the  second-order  terms  are  neglected. 


2.2.2  Pressure  Coefficient 


The  exact  isentropic  pressure  coefficient  can  be  expressed 
as 

2  f  r-1  1  / 1-\  •> 

CP  ~~  im  \[l  ~  -2^2JW)2  +  2Q  .-1)]  -1  }  (2.32) 

oo 

where  t’  is  the  nondiiens ional  time  in  the  wind-fixed  system. 
The  partial  derivative  with  respect  to  t’  as  measured  by  an 
observer  moving  with  the  harmonically  oscillating  body-fixed 
system  is  given  by 

aV  =  at  '  g  ^(t)rcose—  +  g* (t)  (cose^  -  ?  sm©^) 

(2.33) 

Note  that  Eqs.  (2.22),  (2.31)  and  (2.33)  contain  Revell’s 

equations  (Ref.  21)  for  rigid-body  oscillations  as  a  special 
case . 

Substituting  Eqs.  (2.33)  and  (2.28)  into  (2.32)  and,  after 
binomial  expansion,  collecting  terras  of  order  <S0  the  pressure 
coefficient  for  unsteady  flow  C*  can  be  written  as 

=  -  2S0{*1X(1  +  *0x)  -  (1  +  *0x)g’’R#0x 

+  ^Or^ir  +  S’)  +  ik  («^or  +  ~  «'M0x)}  (2.34) 

at  r  =  R ( x ) 

were  S0  is  given  by  Eq.  (2.18)  and  the  pressure  coefficient  for 
the  mean  flow  C°  by  Eq.  (2.16) 


Eqs.  (2.23),  (2.31)  and  (2.34)  together  with  the  Mach  wave 


conditions  at  the  apex,  Eq.  (2.9),  constitute  the  formal 
formulation  according  to  the  body-fixed  coordinate  system. 
However,  to  solve  Eq.  (2.23)  is  rather  tedious.  For  this  reason, 
a  justification  for  simplified  form  of  this  equation  is  sought. 
Van  Dyke  (Ref.  10)  has  shown  that  his  "first-order"  steady  cross 
flow  equation, 


(  1  -  M2)  <fiixx  +  <t>lrr 


-  ^  rf,  =  o 


is  superior  in  yielding  better  results  to  the  linearized 
equation,  which  contains  one  extra  term,  in  the  higher  Mach 
number  range.  For  bodies  in  low-frequency  oscillations,  Platzer 
et  al.  (Ref.  8)  and  Tobak  et  al.  (Ref.  9)  adopted  the  first-order 
equation  formulation  to  obtain  stability  derivatives.  Their 
results  were  found  in  good  agreement  with  computed  results  using 
Euler’s  equations.  Thus,  extending  Van  Dyke’s  first-order 
equation  further  in  the  general  frequency  domain  amounts  to 
neglecting  the  interaction  terms  G1(g,^0)  in  Eq.  (2.23).  In  this 
way,  the  linearized  wave  equation  is  employed  for  the  present 
formulation.  Admittedly,  this  level  of  approximation  is 
i n j us t i f i ab 1 e  mathematically.  Nevertheless,  the  interaction 
effect  due  to  the  term  Gt  can  be  recovered  formally  by  a  Green 
function  approach  as  proposed  in  Chapter  4. 

2.3  Pseudo-Wind-Fixed  Coordinate  System 
The  present  coordinate  system  is  a  hybrid  one.  The  x-axis 
chosen  here  is  not  the  body  axis,  but  remains  rigid  in  motion, 


thus  avoiding  the  complexity  of  using  the  curvilinear  j 

coordinates.  Meanwhile,  the  front-end  of  the  bending  mode  g(x)  | 

| 

when  expressed  in  this  coordinate  is  required  to  attach  to  the 

i 

origin  of  the  x-axis  at  all  times,  as  indicated  in  Fig.  3.  In 

this  way,  the  apex  singularity  can  be  totally  removed.  < 

Admittedly,  such  a  formulation  for  singularity  removal  is  not 

totally  justifiable  but  intended  to  serve  as  a  regular 

approximation.  | 

Hence,  the  linearized  equation  can  be  obtained  in  the  same 
form  as  Eq.  (2.23)  by  letting  g(x)=x-x0  in  Eq.  (2.24),  and  is 
replaced  by  G2(g,^0)  as  follows  < 

G2(g,*0)  =  M2  {2*0xr  +  k2  [  r  <tQx  -  (x-x0)  *0r  +  2ik 

tor  -  r*Oxx  +  (x_xo)'#Orx)  (2.35) 

i 

where  x0  is  a  point  chosen  on  the  x-axis  such  that  p(0)=0.  Since 
this  equation  is  a  degenerate  form  of  the  body-fixed  equation, 
the  justification  for  neglecting  G2  terms  also  remains  valid 
here.  The  mode  shape  g(x)  expressed  in  this  system  now  becomes 
p ( x ) ,  where 

i 

p(x)  =  g(x)  -  ( x  q  —  x  )  (2.36) 

The  potential  is  split  into  two  components,  1 

<t> !  ( x , r )  =  *(x,r)  +  \(x,r)  (2.37) 


i 


17 


where  <t  and  X  represent  the  rigid  part  (x-axis  motion)  and  the 
flexible  part  (p(x)  motion)  of  the  potential,  respectively. 

Thus,  the  boundary  conditions,  after  Tay  1  or-ex pans i on  transferal 
to  the  mean  surface,  become,  at  r=R(x) 

*r  "  R’  =  -  l-ik(x-x0  +  RR’)  at  r=R(x)  (2.38) 

xr  -  R’  xx  =  -  p’(x)  (1  +  0Ox  -  ikp(x) 

+  p(x)0orr  -  p(x)R’0Oxr  at  r=R ( x )  (2.39) 

Notice  that  Eq.  (2.39)  is  essentially  the  same  as  Eq.  (2.14) 
of  the  wind-fixed  system  with  g(x)  replaced  by  p(x).  The 
pressure  coefficients  C£r  and  C*f,  corresponding  respectively 
to  v  and  X,  read 

Cpr=  -  2S0 • { [ ( 1  +  *0X)*X  +  +  D 

+  ik(*  -  Rtf0x  +  (x-x0)^0r)]}  at  r  =  R(x)  (2.40) 


and 


cy  =  -  2S0  •  { ( 1  +  *0x)Xx  +  *0rxr  +  ikv  “  P(*) 

[*Oxr  +  *Ox*Oxr  +  ^or^Orrl)  at  r=R(x)  (2.41) 

where 

C1  =  C1  r  +  C1  f  (2.42) 

p  p  p 
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In  the  present  system  p(0)=0.  Therefore,  this  condition 
guarantees,  all  second-order  coupling  terms  in  the  bracket  of  Eq. 
(2.41)  when  multiplied  by  p(x)  become  necessarily  finite  at  the 


apex . 

It  should  be  noted  that  although  the  apex  singularity  is 
removed  in  the  pseudo-wing-fixed  system  Eq.  (2.41)  is  still  not 
uniformly  valid,  because  these  second-order  derivative  terms  will 
bear  local  singularities  at  other  body  slope  discontinuities  such 
as  the  cone-cylinder  junctions.  In  the  case  of  the  body-fixed 
system,  however,  such  singularities  do  not  appear,  since  the 
unsteady  pressures  only  involve  first-order  derivative  terms. 


2.4  Generalized  Forces  and  Stability  Derivatives 
One  of  the  main  purposes  of  the  unsteady  aerodynamic 
computations  is  to  determine  the  generalized  forces  (Q,j),  which 
appear  in  the  Lagrange’s  equations  of  motion  for  the  structural 
system.  There,  the  generalized  force  Qw  represents  the  work 
done  by  the  aerodynamic  force  FJ  per  unit  displacement  g1  with 
all  other  generalized  coordinates  fixed. 

In  the  present  case,  once  the  unsteady  pressure  coefficients 
as  provided  by  Eqs.  (2.17),  (2.34)  or  (2.42)  are  determined,  the 

generalized  forces  can  be  computed  according  to 


Q 


i  j 


-1 

S  r  e  f 


,1 


0  0 


C-  <  }  > 

P 


R(g<  1  > 


+  RR ’ g ’ ( 1 > ] cos2 ©d©dx 

(2.43) 


where  g '  1  >  is  the  Ith  structural  mode;  C £ ( J ■  is  the  pressure 
coefficient  due  to  the  Jth  mode  of  motion  and  Sr e f  is  the 


reference  area  of  the  body.  Here,  Sref  is  defined  as  the  based 
area  for  the  open-end  bodies  and  the  maximum  cross-section  area 
for  the  closed-end  bodies. 

As  will  be  shown  in  Chapter  5,  the  generalized  forces 
resulting  from  the  unsteady  flow  caused  by  the  changing 
deformations  of  the  body  surface  occur  in  the  equations  of 
flutter.  Thus  the  determination  of  reliable  flutter  boundaries 
depends  on  the  accurate  evaluation  of  the  generalized  forces. 

Once  the  generalized  forces  have  been  determined  the 
stability  derivatives  for  rigid  low-frequency  body  motion  can  be 
determined  in  the  following  form 


CNi  =  -RE<Ql2) 


CMi  =  -RE(Q22) 


C  Na  +  C  N  q  =  -IM(Q12  )/k 


CMa  +  CMq  =  "IM(Q22>/k 


C. .  =  -RE(Q, , )  /k2 
La  11 

cMi  =  -RE(Q21 )/k2 


(2.44) 


where  ”1"  is  the  plunging  mode  and  ”2"  the  pitching  mode,  and  RE 
means  the  real  part  and  I M  the  imaginary  part. 
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CHAPTER  3 


STEADY  MEAN  FLOW 

The  purpose  of  solving  the  steady  mean  flow  problem  is  two¬ 
fold.  First,  it  is  necessary  to  obtain  the  mean  flow  velocities 
on  the  body  in  order  to  proceed  with  the  unsteady  flow 
computations  because  the  unsteady  flow  potential  and  the  unsteady 
pressure  coefficient  depend  on  the  mean  flow  solution.  Second, 
it  is  important  to  establish  a  robust  unsteady  computational 
procedure,  which  can  be  first  established  by  the  steady  mean  flow 
studies  as  a  primary  step,  as  the  steady  mean-flow  is  provided 
with  more  readily  available  data  for  result  verifications. 

It  has  been  noted  that  first-order  theory  predicts  the  cross 
flow  more  accurately  than  the  axial  flow  because  smaller 
disturbances,  in  the  convective  direction,  are  involved.  This 
suggests  that  it  is  more  important  to  refine  the  axial  flow  than 
the  cross  flow.  In  fact  this  type  of  hybrid  theory  has  been 
suggested  previously  by  Van  Dyke  (Ref.  10)  involving  a  linear- 
equation  model  for  the  cross  flow  and  a  nonlinear-equation  model 
for  the  mean  flow. 

For  the  cross-flow  model  the  use  of  this  hybrid  theory  can 
only  be  justified  by  previous  numerical  studies.  In  the  past  Van 
Dyke  (Ref.  10)  obtained  accurate  results  for  bodies  at  steady 
angle  of  attack,  and  Tobak  and  Wehrend  (Ref.  9)  obtained  results 
for  the  static  and  dynamic  stability  derivatives  for  cones,  which 
are  in  good  agreement  with  those  obtained  based  on  the  Euler’s 
solution.  Thus,  in  this  chapter,  together  with  the  solution  of 
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< 


the  linear  equation  for  the  mean  flow,  Eq.  (2.7),  a  nonlinear 
equation  for  the  mean  flow  is  derived.  Its  solution,  by  I 

employing  an  iterative  procedure  first  suggested  by  Van  Dyke,  is 
presented.  In  the  latter  development  and  for  legends  in  Figures 
presented,  this  unsteady  hybrid  model  is  termed  "nonlinear”.  ( 

3.1  Solution  of  the  Linear  Equation 

The  governing  equation  for  the  mean  flow  as  derived  in 
Chapter  2,  Eq.  (2.7)  can  be  written  as: 

(1  *  M2)*0xx  +  *0rr  +  i-  *0r  -  0  (3.1) 

Along  the  apex  Mach  wave  and  on  the  body  surface  the 
boundary  conditions  are: 

'  *or  =  0  at  x  -  pr  <_  0  (3.2) 

<t>  0r  “  H’  0Ox  =  at  r  =  R  ( x )  (3.3) 

Since  Eq.  (3.1)  is  linear,  it  admits  superposition  of 
solutions.  Therefore,  for  a  body  of  revolution,  whose  upstream 
pointed  end  is  at  x  =  0,  the  solution  to  Eq.  (3.1)  can  be  sought  as 
a  superposition  of  supersonic  sources  along  the  x-axis.  The 
potential  at  (x,r)  can  be  expressed  as  (Ref.  13) 

p*-£r  f ( f ) 

*0(x,r)  -  -  -  t{-  -  df  (3.4) 

J0  /(  x-  I )  2  r2 

where  i  denotes  the  distributed  source  location  along  the  x-axis 


and  f(<)  is  the  source  strength  at  t. 


The  velocity  components,  after  differentiation  of  Eq.  (3.4), 
are  given  by 

/3r  f  ’  (  t )  d; 

/(  x- 1 )  2  ~  p2  r2 

x~fir  (x-Qf’(f)d  l 

/(x-l  )2-p2r2  (3.5) 

Eqs.  (3.4)  and  (3.5)  satisfy  the  conditions  given  by  Eq. 
(3.2)  if  f ( 0 ) =  0  and  the  source  strength  function  f  will  be 
determined  from  the  tangency  condition,  Eq.  (3.3).  This  leads  to 
a  Volterra  integral  equation,  which  can,  in  general,  only  be 
solved  numerically.  The  procedure  introduced  by  von  Karman  and 
Moore  (Ref.  13)  proceeds  by  replacing  the  body  by  one  consisting 
of  a  head  cone  and  a  sequence  of  truncated  cones.  The  procedure 
is  clearly  described  in  the  standard  texts  such  as  Sauer  (Ref. 

22)  and  Ferri  (Ref.  23).  Its  generalization  to  the  unsteady 
problems  is  described  in  the  next  chapter. 

3.2  The  Nonlinear  Equation  and  Its  Solution 
A  nonlinear  equation  for  the  mean  flow  perturbation 
potential  ^0(x,r)  can  be  obtained  by  starting  with  the  full 
potential  as  Q(x,r,©)=x+0o(x,r).  Substituting  this  expression 
into  Eqs.  (2.2)  and  (2.3)  yields 


*0r  (x,r)  =  p  | 

0 


*0x  r) 


»> 
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+*  <pn  r  r  +  -  0,~.  “  M2  f* - ^ v  v  $r  r-  **“  —  0n  _  ) 

oo  uxa  '  «j  r  r  r  go  l  '  u  x  x  ’  u  r  r  j^  ur 

(2«*0X  +  *0x  +  *0r  ) 

+  2<*0r  ^Oxr  (  1  +  ^Ox  -*  (3.6) 

-t-oj  (  2  $  .  +  oj  2  )  +  ^  2  ^  T 

^Oxx'^Ox  *^0s  '  "Or'OrrJ 

According  to  Van  Dyke  (Ref.  10),  Eq.  (3.6)  can  be  solved  by 
an  iterative  procedure.  Letting  L(*)  be  the  left-hand  side 
operator  and  R(*)  the  right-hand  side  operator,  Eq.  (3.6)  can  be 
simply  expressed  as 


L(*0)  =  R(*0) 


(3.7) 


where  the  potential  ~$Q  is  the  first-order  iteration  obtained  by 
solving  the  linearized  equation  L(?o)=0.  Making  use  of  this 
fact,  the  right-hand  side  R  (  )  can  be  expressed  as 

R(^0)  =  M«(  ^oxr^or  +  [2  +  (r  -  DM*]#OXX#OX  +  ^0rr^2r(3.8) 


r-1  -  -  r-1 

—  — g- — M^  ^oxx^Or  +  ^Onr^Ox^Or  +  (  1  +  — g — )  ^Oxx^Ox^ 


For  the  second  order  iteration  all  the  triple  products  in 
Eq.  (3.8)  can  be  neglected  except  the  term  ?0rr^ri  which, 
as  can  be  shown  from  the  slender-body  theory,  is  of  the  same 
order  as  the  other  two  quadratic  terras.  Thus,  the  equation  for 
the  second-order  iteration  can  be  written  as 
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L^0  >  -  Mi  l^ox^or  +  [2  +  ^  0  x  x  ^Ox  +  ^orr^orl 

(3.9) 


To  solve  Eq.  (3.9)  Van  Dyke’s  particular  solution  <»0  is  employed. 
*0  =  MLC#0x(#0  +  Nr  *Or>  -  J  ^0  r  1  (3.10) 


where  N=  (  f+ 1 ) M^/ [ 2 (M^-l ) ] .  It  should  be  noted  that  , 

when  substituted  on  the  left-hand  side  of  Eq.  (3.9),  besides 
accounting  for  all  the  terms  on  the  right-hand  side,  also  gives 
some  triple  products  involving  x  derivatives  of  f Q  .  These 
triple  product  terms  are  of  order  equal  to  or  higher  than  those 
terms  already  neglected  in  Eq.  (3.8).  Hence,  for  the  second 
order  iteration,  v0  can  be  considered  as  an  exact  particular 
solution.  Then,  <fiQ  can  be  expressed  as  0O  =  XQ+(»O,  where  XQ 
satisfies  L(Xo)-0. 

Both  potential  <fiQ  and  ~$Q  must  satisfy  the  Mach  wave 
condition  Eq.  (3.2),  this  dictates  that  both  / 0  and  XQ  should 
also  satisfy  the  same  condition  independently.  The  problem  of  X0 
is  just  the  first  order  problem,  with  the  tangency  condition  Eq. 
(3.3)  replaced  b y 
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< 

X0r(x,r)  -  R’X0x(x,r)  =  R  ’  (  1  +  <*0  x  (  x  ,  r  )  )  (3.11) 

-  K0r(x,r)  at  r  =  R ( x )  ( 

The  complete  second  order  perturbation  velocities  are  found 
as  the  sums  of  the  contributions  from  *0  and  X0 .  Then,  the 
pressure  coefficient  can  be  calculated  according  to  Eq.  (2.16). 

3.3  Results  and  Discussions 

To  demonstrate  the  present  method  and  to  validate  its  1 

procedures  several  computed  cases  are  presented:  the  mean  flow 

total  velocity  for  a  cone-cylinder,  mean  flow  pressures  for  a 
parabolic-ogive,  and  an  ogive-cy 1 inder-boattai 1  bodies  are  shown  l 

as  examples. 

Figure  4  demonstrates  that  the  present  linear  method  yields 
correct  values  for  the  velocity.  It  is  seen  that  the  present 
result  is  in  good  agreement  with  the  USSAERO  result  as  well  as 
those  obtained  by  other  theories  (see  Ref.  10)  for  a  cone- 
cylinder  body.  The  deviation  of  USSAERO  result  on  the  aft- 
cylinder  is  probably  caused  by  an  erroneous  wave  influence 
generated  by  the  lower  junction  of  the  cone-cylinder.  In  Figs.  5 
and  6  the  linear  and  the  nonlinear  results  are  compared  with 
those  computed  by  USSAERO  code,  and  by  the  exact  method  of 
characteristics  (Ref.  24),  for  a  26%  thick  ogive-cylinder  body  at 
Mach  numbers  of  2.0  and  3.0,  respectively.  The  nonlinear  results 
compare  very  well  with  those  computed  by  the  exact  method  of 
characteristics.  It  is  seen  that  the  nonlinear  effect  due  to  the 
thickness  is  substantial  from  the  apex  to  raid-body. 
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Next,  the  hybrid  theory  (nonlinear  iterative  scheme  for  the 
mean  flow  and  linear  for  the  cross  flow)  is  applied  to  a  16% 
tnick  og 1 ve-cy 1 i nde r-b oa 1 1 a i  1  body  at  Mach  number  of  3.0  and 
placed  at  moderate  angles  of  attack  (aQ=3.2  and  6.3).  Again, 
very  good  correlations  are  found  with  the  computed  results  of  the 
Parabolized  Na v i er-S t okes  (PNS)  code  and  the  Euler  code  (Ref.  25) 
for  both  cases  in  Figs.  7a  and  7b.  Considerable  deviations 
between  the  linear  and  the  nonlinear  results  are  again  observed 
particularly  on  the  windward  side  of  the  ogive  part  of  the  body. 

It  can  be  concluded  that,  as  long  as  the  flow  remains 
attached,  the  present  nonlinear  method  should  yield  results  in 
favorable  agreement  with  those  obtained  by  computational  methods 
in  the  supersonic  Mach  number  range.  This  agreement  also  implies 
that  for  the  given  range  of  Mach  number,  body  thickness  and 
angle-of-attack ,  effects  of  r o t a t i ona 1 i t y  as  introduced  by 
supersonic  shock  waves  are  nearly  inconsequential. 
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CHAPTER  4 


UNSTEAD i  FLOW  COMPUTATIONS 

In  this  chapter,  the  method  of  solution  for  a  uniform  flow 
over  a  body  performing  bending  oscillations  is  presented.  The 
formulations  in  the  wind-fixed,  body-fixed  and  pseudo-wind-fixed 
coordinate  systems,  as  described  in  Chapter  2,  are  used.  To 
compute  the  oscillatory  flow  in  the  cross  plane,  a  line  doublet 
distribution  scheme  (see  Fig.  8)  is  adopted.  The  strength  of  the 
doublet  along  the  axis  is  modeled  according  to  the  Harmonic 
Gradient  model  for  treating  wing  planforms  (Ref.  4)  This  model 
is  capable  of  rendering  the  unsteady  potential  solution  and  its 
convective  gradient  uniformly  valid  throughout  the  complete 
frequency  domain. 

To  verify  the  present  method,  numerical  examples  for  various 
body  shapes  are  presented,  in  terms  of  unsteady  pressures, 
stability  derivatives,  generalized  forces,  and  aerodynamic 
damping,  and  compared  with  various  theories  and  measured  data. 

4.1  The  Integral  Solution 

The  general  integral  solution  to  the  unsteady  wave  equation, 
e.g.  Eq.  (2.8)  or  Eq.  (2.23),  can  be  obtained  by  applying  the 
Green  function  method. 

_i  pX-Ar  a 

4>i  ( x  ,  r )  =  ^  j  F(  f  )yp  K(x-«,/Jr)di  +  *g  (4.1) 

o 


where 
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<t>a  =  J  J  Gn  (  (  f  >  7  )  >  2)  K(x-f,  /Sr)d7d« 

A 

The  function  Gn(0o,£)  is  the  inhomogeneous  term  of  the  unsteady 
wave  equation,  Gr  in  Eq.  (2.23)  or  G2  in  Eq.  (2.35)  and  the  area 
A  is  defined  in  the  domain  downstream  of  the  Mach  wave,  emanating 
from  the  apex,  excluding  the  body,  i.e. 

R ( x )  <  7  <  i/p  ,  for  0  <  i  <  1 

0  £  7  <  i/p  ,  for  1  <  I  i  «• 

The  first  term  of  Eq.  (4.1)  represents  the  linear  unsteady 
solution  as  modeled  by  the  distribution  of  doublets  along  the  x- 
axis,  and  the  second  term,  <fig,  the  mean  flow-mode  shape 
interact  ion . 

In  the  present  analysis,  we  shall  drop  the  term  in  Eq. 
(4.1)  for  simplicity.  But  in  principle  <fit  can  be  included  in  the 
analysis  since  <j>0  is  known  from  the  mean  flow  computation  and 
g(x)  is  given.  The  kernel  function  K  is  an  elementary  solution 
of  Eq.  (2.6)  (see  Garrick,  Ref.  26) 

K(x-f,  pr )  =  e-i/t(x~f)  cosvR  (4.2) 


where  R  is  now  the  hyperbolic  distance  and  is  defined  as 


R  =  /( x-  f  )2  -  p2  r2  ,  ft  =  kM \jp*  and  V  =  kM  J  p2 


and  F(f)  is  the  dipole  strength  to  be  sought. 


Now,  integrating  Eq.  (4.1)  by  parts,  and  making  use  of  the 


Mach  wave  condition  at  x=^r,  yields 


.  x  ,  r 


-1 
2  it 


dx, 


[f ( x-x0  )  e  1/‘X°]|?S(x0  ,/8r)dx0  (4.3) 


where 


l?S(x0,/rr)  =  §p  I 


x „  cosXR 


pr 


<0  =  x-<  and  =  <Jz2-/32r2 


(4.4) 


Notice  that  from  this  point  onward  "x0"  is  used  to  represent 
the  relative  dipole  coordinate. 

4.2  Harmonic  Gradient  Model 

The  solution  of  Eq.  (4.3)  is  based  on  a  line-doublet  panel 
method  similar  to  Karman  and  Moore’s  (Ref.  13)  procedure  for  the 
source  solution  formulation  to  the  mean  flow.  Thus,  a  set  of  N 
points  with  coordinates  (xjfrj(0)  j=l,...N  are  distributed  on 
the  body  surface,  such  that  Xj , 1 >x^ .  These  points  are  called  the 
control  points.  To  determine  the  induced  potential  at  each 
control  point  the  intersection  of  the  inverse  Mach  cone,  from  the 
control  point  (Xj.rj.O),  with  the  body  axis  is  first  determined, 
(see  Fig.  8).  The  set  of  points  so  obtained  are  given  by 
tj*i=xj~Prj  j  =  l,...N,  with  being  the  body  apex.  The  segment 
between  each  two  points  ( i.  ,  + : )  is  called  a  panel.  Each  panel 
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is  assumed  to  have  different  dipole  strength  Fj  (i).  The 
potential  <j>x  at  (Xj  ,  r_j  )  can  then  be  expressed  as 


-1  f  fxj  +  1 
•rJ>=  27  [  j 


C^i  (xj  -x0  )e_i^X°] 


i  =  1  ~  j  i 


(x0>^rj  )dx0  # 


(4.5) 


In  order  to  achieve  computational  accuracy  and  effectiveness 
for  handling  solutions  in  the  high-frequency  range,  it  is 
important  to  render  the  doublet  solution  and  its  convective 
gradient  uniformly  valid  throughout  the  complete  frequency 
domain.  This  is  to  say  that  the  characteristics  of  the  doublet 
solution  should  be  spatially  harmonic.  Therefore,  the  panel  size 
from  fl  to  iiM  is  regulated  and  maintained  compatible  to  the 
wave  number  generated  along  the  body  in  oscillation  (see  Fig.  9). 
This  is  the  Harmonic  Gradient  concept  introduced  by  Chen  and  Liu4 
for  unsteady  supersonic  computations.  Following  this  concept, 
one  can  model  the  integrand  of  Eq.  (4.5)  in  a  similar  manner, 
i  .  e  .  , 

|^-  [Fi  (xJ-x0  )e~1/iX°]  =  [ai(Xj-x0)  +  bje  1/1X°  (4.6) 


where  at  and  bt  on  the  right-hand  side  are  complex  constants 


4.3  Evaluation  of  Velocities 


With  the  Harmonic  Gradient  (H-G)  model  of  Eq.  (4.6),  Eq. 
(4.5)  can  be  differentiated  to  obtain  the  discretized  velocities 
x  and  ^  r  ,  i  .  e  . 


* 


(x 


T;  > 


e_i^xo 


)dx0 


(4.7) 


<t> ,  r  (  x  j  ,  r  j  )  = 


-1 

(  1 

r*j-f 

2  IT 

l  j 

x  —  ? 

1 

=  1 

J 

as 

dr? 

S  (x0  ,  /Jr 

[ a i ( x j -x0  )  +  bt  ]e 


-i  MX, 


(4.8) 


For  detailed  evaluation  of  these  integrals  one  is  referred 
to  Appendix  B.  There  it  is  shown  that  Eqs.  (4.5),  (4.7)  and 

(4.8)  can  be  integrated  numerically  by  employing  Lashka’s  (Ref. 
27)  exponential  substitution  scheme. 

4.4  Method  of  Solution 

After  Eqs  (4.5),  (4.7)  and  (4.8)  have  been  integrated,  the 

potential  and  the  velocities  at  (x; , r ^  )  can  be  expressed  as 
functions  of  the  as  ’s  and  bs  ’s.  However,  the  b4  ’s  can  be 
determined  as  functions  of  the  ai ’s  by  imposing  the  condition 
that  the  potential  must  be  continuous  between  each  two  adjacent 


panels.  With  this  condition  the  following  recurrence  formula  can 


be  established,  i.e. 

bn  =  -T 

n  i  v  L  i/i 

j  =  2 

+  £a-(  1  -  e1^in-i  )  n  >  1 

i  V  ’ 

b:  -  0 


(4.9) 


where  b^O  is  due  to  the  application  of  the  apex  Mach  wave 
condi t ion . 

Since  the  bs ’s  are  expressed  in  terms  of  a^s,  the  latter 
can  be  evaluated  by  applying  the  tangency  condition  at  the 
control  points.  In  matrix  form  this  condition  can  be  expressed 
as 


[W3i]  {a*}  =  { B  j }  (4.10) 

where  [Wj;]  can  be  expressed  as  [ v i  t  ] -H * ( x i ) [ U j . ]  and  Vji  and  Uj ; 
are  the  velocity  influence  coefficients  in  r  and  x  directions 
respectively  of  the  oscillatory  flow  at  x^  ,  r ^  due  to  the  panel 
i.  The  r i gh t -hand-s i de  of  the  tangency  condition  evaluated  at 
the  control  point  j  is  denoted  by  { B ^ } ,  which  represents  the 
given  downwash. 

In  supersonic  flow,  the  governing  equation  is  hyperbolic, 
the  problem  becomes  an  initial  value  problem;  hence  the  matrix 
[W  , ]  is  a  lower  tridiagonal  one.  The  solution  method  for 


solving  the  system  of  equations  given  by  Eq.  (4.10)  is 
straightforward. 

Once  the  { a t }  are  determined,  we  can  compute  the  velocities 
and  the  potential  at  the  control  points  and  therefore  the 
unsteady  pressure  coefficient.  Once  the  unsteady  pressures  have 
been  determined,  the  generalized  forces,  can  be  determined  by 
using  Eq.  (2.43). 


4.5  Results  and  Discussions 
To  verify  the  present  method,  numerical  examples  are 
presented  in  terms  of  unsteady  pressures,  stability  derivatives, 
and  generalized  forces.  Free-free  mode  aerodynamic  damping  for 
bodies  in  bending  oscillations  are  presented  for  various 
configurations,  including  that  of  the  Saturn  SA-1  launch  vehicle. 

4.5.1  Results  According  to  Various  Coordinate  Systems 

Figures  10,  11  and  12  present  the  in-phase  and  out-of-phase 

pressure  coefficients  of  a  10S>  thick  cone  at  Mach  number  M  =2.0 

oo 

and  reduced  frequency  k=2.0  in  the  wind-fixed,  body-fixed  and 
pseudo-wind-fixed  coordinates.  The  oscillating  cone  performs  in 
pitching  mode,  f i r s t -b end i n g  mode  and  second-bending  mode, 
respectively.  Free-free  mode  beam  theory  was  used  to  determine 
these  modes.  For  the  rigid  mode  oscillation  in  Fig.  10  the 
pseudo-wind-fixed  and  the  body-fixed  systems  become  identical; 
thus  only  one  result  is  presented. 

In  general,  the  results  of  the  body-fixed  and  the  pseudo¬ 
wind-fixed  systems  and  those  of  the  slender  body  theory  (Ref.  6) 
are  in  good  agreement.  In  contrast  to  these  results,  the  in- 


34 


phase  pressures  of  the  wind-fixed  results  persistently  show  the 
effects  of  the  apex  singularity  in  all  cases,  as  expected. 
Consequently,  the  overall  pressure  distributions  downstream  are 
contaminated  by  this  singular  behavior  originated  from  the  apex. 
It  is  noted  that  when  the  oscillation  center  x5  is  placed  at  the 
apex,  or  the  mode  shape  g(x)  at  the  apex  is  zero,  the  apex 
singularity  disappears  and  all  wind-fixed  results  are  in  close 
agreement  with  the  others.  Also,  it  can  be  observed  that  the 
pressures  due  to  the  flexible  modes  are  one  order  higher  than 
those  of  the  rigid  mode,  and  that  the  out-of-phase  pressures 
resemble  the  mode  shape.  The  reason  for  this  latter  behavior  can 
be  simply  analyzed  from  the  slender  body  limit.  There,  the  out- 
of-phase  pressure,  as  shown  in  Appendix  C,  is  proportional  to 
(R(x)g(x))’  and  because  R(x)  goes  to  zero  at  the  apex,  the 
dominant  term  is  R’(x)g(x). 

From  Figs.  13  through  18,  generalized  aerodynamic  forces  on 
a  cone,  a  parabolic-ogive  and  a  cone-cylinder  oscillating  in 
first  and  second  bending  modes  with  reduced  frequency  of  k=2.0 
are  plotted  versus  freestreara  Mach  number.  (Due  to  the 
singularity  originated  from  the  apex,  the  generalized  forces 
resulted  from  the  wind-fixed  coordinate  system  will  not  be 
presented  here.)  As  expected,  the  present  results  approach  the 
slender  body  results  in  the  decreasing  order  of  thickness. 

While  the  pseudo-wind-fixed  and  body-fixed  results  are  in 
satisfactory  agreement,  their  deviations  increase  with  increased 
Mach  number  and  thickness.  In  Fig.  18,  the  sudden  departure  of 
the  pseudo-wind-fixed  phase  angle  (argument)  could  be  caused  by 


35 


the  mode-shape /'expansion-fan  interaction  which  is  further 
amplified  by  the  second  derivative  term  <j>  0  r  r  .  In  general,  it  can 
be  observed  that  again  higher  order  modes  result  in  higher  value 
of  generalized  forces.  For  example,  the  magnitude  (modulus)  of 
the  generalized  forces  for  the  second  bending  mode  is  about  twice 
that  of  the  first  one. 

In  terms  of  the  effects  of  Mach  number,  it  can  be  seen  that 
when  the  Mach  number  approaches  the  low  supersonic  regime,  i.e. 

M  =1.5,  all  magnitudes  increase  rapidly.  This 
trend  is  similar  to  that  obtained  for  rigid  modes.  When 
approaching  the  higher  Mach  number  range,  the  changing  rate  of 
the  force  magnitudes  appears  to  be  less  sensitive  to  the  Mach 
number  for  a  given  body  thickness  or  given  mode  shape.  Similar 
trends,  in  the  high  supersonic  Mach  number  range,  are  found  in 
the  solution  to  Euler’s  equations  for  steady  flow. 

To  simplify  the  matter,  in  the  following  sub-sections,  only 
results  in  tha  body-fixed  coordinate  system,  are  presented. 

4.5.2  Effects  of  Frequency 

In  the  low  -frequency  limit,  the  damping-in-pitch  moment 
coefficients  for  a  parabolic-ogive,  an  ogive-cylinder  and  a  cone- 
frustrum  body  are  presented  in  Figs.  19,  20  and  21.  Throughout 
the  supersonic  range  the  present  results  are  found  in  good 
agreement  with  those  of  Platzer’s  (Ref.  8)  linearized  method  of 
characteristics  ( LMOC )  ,  Tobak  and  Wehrend’s  (Ref.  9)  cone  theory, 
which  are  limited  to  the  low  frequency  domain,  and  various 
experimental  data  (Ref.  8).  Across  the  frequency  range,  the 
magnitude  and  the  phase  angle  of  the  generalized  forces  Qrj  for  a 
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cone  and  an  ogive  in  pitching  motion  are  presented  in  Figs.  22a 
and  22b.  The  force  magnitude  increases  with  increased  frequency, 
whereas  the  phase  angle  is  less  frequency  dependent.  Notice  also 
that  the  effects  of  hody  shape  become  apparent  in  the  high 
frequency  range.  All  results  merge  in  the  low  frequency  limit. 
Similar  trends  were  found  in  the  sonic-flow  studies  by  Landahl 
(Ref.  28). 

In  the  high  frequency  limit,  Fig.  23  presents  results  of  the 
present  method  compared  against  the  p i s t on- t heo ry  (Ref.  29) 
results  for  a  very  slender  parabolic  ogive.  The  thickness  ratio 
r=.02  for  this  case  is  selected  baased  on  the  order  analysis 
tM<>)>k«l  and  M^kil.O  as  required  by  the  piston  theory.  For  M  =  1.5 
the  agreement  seems  to  be  very  good  for  the  two  selected  reduced 
frequencies,  k  =  4.0  and  7.5.  It  is  also  interesting  to  compare 
the  effects  of  frequency  and  flow  dimensionality  on  unsteady 
pressures.  Figure  24  shows  comparisons  of  unsteady  pressure 
coefficients  for  a  5.7*  cone  and  a  flat  plate  pitching  about 
the  apex  at  M  =2.0.  The  flat  plate  results  are  computed  by  the 
LPP  code  (Ref.  30).  As  expected,  the  unsteady  pressure  magnitude 
for  an  oscillating  cone  is  smaller  than  that  of  a  flat  plate  at 
k=1.0  and  k=2.0.  Similar  to  the  case  of  steady  supersonic  cone 
and  wedge  flow,  the  present  finding  shows  that  the  cone  in 
oscillation  yields  weaker  compression  than  the  flat  plate  as  a 
result  of  the  three  dimensionality  of  the  flow,  irrespective  of 
the  oscillation  frequency. 
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4.5.3  Nonlinear  Rerults 


The  effects  of  the  nonlinear  mean-flow  in  the  unsteady- 
aerodynamic  forces,  are  presented  in  Figs.  27  and  28.  In  these 
figures  the  legend  ’’Present  HPP”  represents  the  HPP  linear 
results  to  distinguish  them  from  the  "Present  HPP  nonlinear" 
results.  It  should  be  noted  that  the  latter  results  are  obtained 
based  on  the  unsteady  hybrid  approach  developed  in  Chapter  3. 

Figure  27  shows  that  a  good  agreement  of  HPP  nonlinear 
results  with  Brong’s  (Ref.  31)  exac  ;  Euler  unsteady  results  is 
obtained  for  the  pitch-damping  forces  for  a  cone.  The  linear 
results  deviate  from  the  nonlinear  ones  as  the  Mach  number 
increases.  The  nonlinear  effect  is  enhanced  by  either  increasing 
the  hypersonic  parameter  M  r  or  by  decreasing  the  Mach  number 
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toward  the  transonic  range. 

Figure  28  shows  that  the  present  HPP  methods  are  in  fair 
agreement  with  the  measured  damp i ng- i n-p i t ch  moment  coefficient 
for  a  20%  thick  ogive-cylinder  throughout  the  Mach  number  range. 
The  computed  results  of  SPINNER  code  (Ref.  32)  and  Ericsson’s 
(Ref.  32)  show  large  discrepancies  with  the  measured  data;  in 
fact,  weak  dependency  on  the  Mach  number  range  was  found  in  these 
results.  By  contrast,  strong  Mach  number  dependency  is  shown  in 
the  results  of  the  HPP  code,  which  show  a  favorable  trend  with 
the  measured  data.  However,  no  appreciable  difference  is  found 
between  the  HPP  linear  and  nonlinear  results  for  this  case. 

From  these  figures,  it  is  seen  that  the  unsteady 
aerodynamics  can  be  altered  substantially  by  the  mean  flow 
influence  through  the  tangency  condition  and  the  pressure 
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coefficient.  These  improved  results  prompt  further  development 
of  the  nonlinear  method  for  aeroelastic  applications. 

4,5.4  Elastic  Bodies 

Next,  computed  results  of  the  HPP  method  for  an  elastic 
cone-cylinder  body  are  compared  with  the  aerodynamic  damping  data 
measured  by  Hanson  and  Dogget  (Ref.  14)  where  the  aerodynamic 
damping  derivative  is  defined  as 

c;  =  2k/z  (  C  A  /Cc  r  )  =  -  ^4^  }  (4.11) 

and  fi  is  the  mass  ratio.  Physically,  this  coefficient,  , 
represents  the  ratio  of  aerodynamic  damping  to  the  critical 
damping.  The  reduced  frequency  of  the  first  bending  mode  (Fig. 
29a)  lies  in  the  range  between  1.12  to  1.6  corresponding  to 
Mw=3.0  to  1.5.  For  the  second  bending  mode  (Fig.  29b),  it  lies 
between  2.9  and  4.2  for  the  same  Mach  number  range.  It  is  seen 
that  the  present  results  establish  close  trends  with  the  measured 
data.  By  contrast,  all  quasi-steady  theories  yield  inferior 
predictions.  Due  to  an  inconsistent  formulation,  in  their 
boundary  condition,  Bond  and  Packard  theory  (Ref.  11)  results  in 
considerable  discrepancy  with  the  measured  data,  as  can  be  seen 
in  Fig.  29b . 

Aeroelastic  analysis  of  the  Saturn  SA-1  launch  vehicle  are 
presented  in  Figs.  .30  and  31.  Steady  mean  flow  pressure,  and  the 
in-phase  and  out-of-phase  pressures  for  the  vehicle  in  rigid  mode 
are  computed  in  Figs.  30a  to  30c.  The  in-phase  and  out-of-phase 
pressures  practically  follow  the  same  trend  as  that  predicted  by 
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the  slender  body  theory.  Clearly  the  deviation  between  the 
slender  body  results  and  the  present  HPP  results  are  due  to  the 
inclusion  of  the  Mach  number  and  the  thickness  effects  in  the 
latter  approach. 

The  mode  shape  determined  by  experiments  at  NASA-Langley 
(Ref.  15)  is  used  as  input  to  compute  the  aerodynamic  damping 
coefficient.  This  coefficient  is  now  defined  as 

Ch =  -  I  m ( Q3 3  ) / n k 7 ,  where  "3"  denotes  the  first  bending  mode  and 
7  is  a  parameter  involving  the  body  shape  and  the  mode  shape 
(7=1.763,  see  Ref.  15).  The  natural  frequency  for  the  actual 
vehicle  is  2.8  Hertz.  However,  because  of  the  large  reference 
length,  the  reduced  frequency  k  lies  between  1.4  and  2.53  for  a 
Mach  number  range  of  3.0  to  1.2,  respectively.  Therefore,  the 
present  case  of  study  needs  an  accurate  prediction  method  in  the 
high-frequency  range.  Again,  good  agreement  is  found  between  the 
present  results  and  the  measured  data.  To  model  this  complex 
configuration,  less  than  100  panels  were  used  in  the  prescribed 
frequency  range.  Consequently,  only  30  seconds  of  CPU  time  in  an 
IBM  3081  were  needed  to  obtain  all  the  data  reported  in  Figs.  30 
and  31. 
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CHAPTER  5 


FLUTTER  COMPUTATION 

One  of  the  main  applications  of  unsteady  aerodynamics  is  the 
computation  of  flutter.  The  phenomenon  of  flutter  is  a  result  of 
an  interaction  of  aerodynamic  forces  and  airframe  response  in 
such  a  way  that  the  structure  receives  energy  from  the  airstream 
rather  than  giving  it  up  as  damping.  In  this  chapter  the  flutter 
boundaries  for  a  7.5*  cone,  in  pitching  and  plunging  oscillations 
are  compared  with  experiments  and  different  theories  for 
different  Mach  numbers.  The  effects  of  the  nonlinear  mean  flow 
are  studied. 

For  the  same  cone  the  divergence  boundary  as  a  function  of 
the  Mach  number  is  presented. 

5.1  Flutter  Equations 

The  equations  of  motion  for  a  body  in  plunging  and  pitching 
oscillations  as  shown  in  Fig.  30  can  be  shown  to  be 

mh(t)  +  Sa«(t)  +  Khh(t)  =  pJI  [Qll^-  +  Qi2a(t)]sref 

sah(t)  +  Ia«(t)  +  Kaa(t)  =  PooU^L[Q21^iil  +  Q22*(t)]sref 

(5.1) 

where  h(t)  and  a(t)  are  'he  generalized  coordinates  for  plunging 
and  pitching  respectively  and  are  defined  positive  down  for  h(t) 
and  clockwise  for  a(t). 
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For  harmonic  oscillations  h(t)  and  a ( t )  can  be  expressed  as 


h  ( t )  =h  e 1  Ult  and  a  ( t )  =  a  e 1  .  Using  the  conventional  notation  for 


the  flutter  equations  we  let  Kh=muj2,  K  =1  u2  ,  p  =  m/rrp  LR(L)2  and 


r2=I  mL2 ,  where  wh  and  w  are  the  natural  frequencies  of  the  body 

in  plunging  and  pitching  oscillations  respectively,  p  is  the  mass 

ratio,  and  r  is  the  dimensionless  radius  of  gyration  about  the 
a 

pitch  axis. 

Eqs.  (5.1)  can  now  be  expressed  in  a  matrix  form  as  follows 


r 


WK  ,  *  ,  a  „  *  _  1_  Q1  1 


1  -  (=*-)  *  (  — )  or  -1-T 

w  uj  2.  p k2 

Ot 


1  Qxjl 

2  pk2 


1  ^ 

2  p k2 


r2  (  1 — )  -  —  ®2-2~  , 

ro.  '  w2  ;  2  pk2  J 


f  h(t)l 


«(t) 


=  {0} 


(5.2) 


In  Eqs.  (5.2)  for  a  fixed  Mach  number  the  unknowns  are  w,  p^ 

and  U  ,  which  represent  the  flutter  frequency,  the  altitude  at 
which  flutter  will  occur,  and  the  flutter  speed  respectively. 

The  U-g  method  (Ref.  33)  is  used  to  determine  the  flutter 
boundary.  Complex  roots  are  obtained  by  introducing  the 
artificial  structural  damping  factor  g,  and  a  root  of  the 
equation  represents  a  point  on  the  flutter  boundary  if  the 
corresponding  value  of  g  equals  zero.  Thus  Eqs.  (5.2)  are 


replaced  by 


♦ 


1 


(=3x.)  >2 

'w  ;  A  2  /ik2 


l 


2  /ik2 


_  I  Qg  2 


2  Mk2 


r2  (1-X2  ) 

CX 


1 fh(t) 


l 


2  /<k2  J 


a  (  t 


=  {0} 


(5.3) 


where  ^2  =  ( )  (1  +  ig). 

If  flutter  exists,  h(t)  and  a ( t )  do  not  vanish  identically. 
Thus,  for  a  nontrivial  solution  to  exist  requires  the  determinant 
of  the  coefficient  matrix  of  Eq.  (5.3)  be  set  to  zero. 
Mathematically  this  amounts  to  solving  an  eigenvalue  problem. 

Solving  the  determinant  for  X  yields  roots  (X1,X2,...)  from 
which  a  new  frequency  and  damping  value  for  each  mode  are 
obtained  as  follows 

W 

a 

"  Real ( \i ) 
g,  =  ( jjM  Imag(  Xt  ) 


The  velocities  corresponding  to  are  obtained  from  k  as 


u .  L 
~ k  ‘ 


Since  the  generalized  forces  are  functions  of  k  the 
eigenvalue  problem  can  be  solved  for  a  number  of  k’s  to  obtain 


43 


( wi  ,gi  , Ui  ) k  for  each  k  and  plot  a  root  locus  (U-g)  diagram.  The 
zero  crossing  of  a  g  locus  denotes  a  flutter  point. 

5.2  Flutter  Results 

For  a  7.5*  cone  at  Mach  numbers  M^=2.0  and  3.0,  the  flutter 
boundaries  in  terras  of  flutter  speeds  and  flutter  frequencies 
versus  the  wind-off  frequency  ratio  are  presented  in  Figs.  31  and 
32.  Two  rigid  modes,  plunging  and  pitching,  are  investigated. 

As  reported  in  Ref.  16,  the  measured  data  were  obtained  in  order 
to  evaluate  the  sufficiency  and  the  applicability  of  the  then 
existing  unsteady  theories  for  flutter  analysis,  namely  the 
quasi-steady  (Q.S.)  approaches  (Refs.  10  and  13)  and  the 
frequency  expansion  theory  (Ref.  7).  For  this  reason,  the 
reduced  frequency  range  is  confined  to  one  below  0.4  so  that  the 
compared  quasi-steady  theories  can  be  valid.  While  all  methods 
yield  rather  accurate  flutter  frequencies,  it  is  seen  that  the 
quasi-steady  method  fails  to  predict  the  flutter  speed 
consistently  with  the  measured  data.  The  present  method  however 
consistently  slightly  underpredicts  the  flutter  speed,  whereas  as 
expected  the  slender  body  theory  predicts  the  most  conservative 
boundary . 

In  order  to  investigate  the  effects  of  the  nonlinear  mean 

flow,  the  flutter  boundaries  for  the  same  cone  are  now  presented 

in  Figs.  33  and  34  versus  the  Mach  number  for  a  wind-off 

frequency  ratio  of  w./w  =1.8.  From  the  comparison 

n  a 

with  the  measured  data  (Ref.  16)  it  is  seen  that  consistent 
improvement  in  trends  are  obtained  over  the  linear  ones. 

However,  the  predicted  boundaries  become  less  conservative  in  the 
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order  of  slender  body,  the  HPP  linear  and  the  HPP  nonlinear 
results. 

The  phenomenon  of  divergence  is  inherent  to  the  static 
ae r oe  1  as t i c i t y .  Divergence  occurs  when  the  static  moment  created 
by  the  aerodynamic  forces  equals  the  elastic  moment 


L  q 


ref 


ac,. 

d  a 


K 

Ot 


In  Ref.  16,  the  divergence  parameter  is  defined  as 
2Lq7rR(L)2  _  2 

k  ac.  ' 

a  -T— 1 - 


In  Fig.  35  the  divergence  boundary  for  a  7.5*  is  presented 
versus  the  Mach  number.  Interestingly,  the  nonlinear  results 
behave  like  the  slender  body  results  in  this  case,  whereby  the 
former  show  little  dependency  on  Mach  number  up  to  M<>)(  =  5.0. 
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CHAPTER  6 


ASYMMETRIC  BODIES 

In  this  chapter,  the  formulation  for  an  asymmetric  body 
performing  bending  oscillations  is  presented.  The  differential 
equations,  boundary  conditions,  and  pressure  coefficients  in  the 
body-fixed  coordinate  system  are  properly  formulated,  and  the 
solution  of  the  first  order  equation,  using  a  spline  panel 
method,  is  obtained.  Comparisons  of  the  steady  and  unsteady 
pressures,  forces  and  moments  for  conical  bodies  of  various  cross 
sections  are  made  with  other  methods,  whenever  available. 

6.1  Formulation 

As  shown  in  Appendix  A,  Eq.  (A. 7),  the  full  potential  in  the 
body-fixed  system  for  asymmetric  bodies  can  be  written  as 

0(x,y,z,t)  =  x  +  hx ( x , t ) z  +  ♦ ( x , y , z , t )  (6.1) 

The  linearized  equation  for  the  perturbation  potential  * 
obtained  by  substituting  Eq.  (6.1)  into  Eqs.  (2.2)  and  (2.3)  is 
given  in  Appendix  A  by  Eq.  (A. 13).  There  it  is  also  shown  that 
the  governing  equations  for  <fiQ  and  <fil  ,  in  cylindrical 
coordinates,  are  given  by 

(i— M*)*oxx  +  *Orr  +  f  *  0  r  +  *oe©  =  0  (6.2) 
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(  1-M2)*lxx  +  *lrr  +  -p  <Plr 

+  FF  ^ee  -  2iM2„*ix+  M2k2^1  =  G  3  (  g ,  0O  )  (6.3) 


where  G3 ( g , )  represents  the  mean  flow  and  flexible  mode 
interaction  (like  in  Chapter  2  Eq.  (2.23)  for  bodies  of 
revolution)  and  is  given  explicitly  as  a  function  of  g  and  <f>0  in 
Appendix  A,  Eq.  (A.  17). 

6.1.1  Boundary  Conditions 

The  condition  at  the  apex  Mach  wave  for  asymmetric  bodies  is 
given  by 


^or  q  ~  ® 

*1  :  =  *ir  =  *ie  =  0 


at  x  <  £r  (6.4) 


and  the  condition  of  the  velocity  to  be  tangent  to  the  body 
surface  at  all  the  times  can  be  expressed  (the  same  as  Eq.  (2.25) 
for  bodies  of  revolution)  as 

(7Q-V8 ) • n  =  0  at  S  =  0  (6.5) 

where  V8  is  the  velocity  of  the  body  surface  and  is  given  by 
V8  =  <5(t)[g’R  cos  ©  ex  -  g  ( cos  ©  er  -  sin  ©  e  )]  (6.6) 
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S  the  body  surface  and  n  the  outward  normal  to  S.  In  the  body- 


fixed  coordinate  system  the  body  surface  S  is  simply  given  by 
S  =  r  -  R(x,e)  =0  (6.7) 

The  normal  n  is  then  obtained  as 

Re  - 

n  =  7S  =  -  Rx  ex  +  er  -  e0  (6.8) 

The  total  velocity  70  is  given  in  terras  of  the  perturbation 
potential  $  in  Appendix  A,  Eq.  (A. 8)  as 

70  =  (1  +  *x  -  hxxz*x)  ex  +  ♦y  ey  +  (hx  +  *2 )  e2  (6.9) 

Now,  expressing  Eq.  (6.9)  in  cylindrical  coordinates  and 
letting 

♦  (  x  ,  r  ,  e ,  t )  =  ( x  ,  r  ,  9)  +  5  ( t )  ^(x.r.e)  (6.10) 

Eq.  (6.9),  when  evaluated  at  the  body  surface,  becomes 

70  =  [1  +  #0x(l  -  5(t)g’ ’Rcose)  +  *(t)#lxlex 

+  [^or  +  5(t)^lr  +  5 ( t ) g ’ cos©] er  (6.11) 

+  t^oe  +  ~  <5  ( t )  g  ’  sine]  eQ  at  r  =  R(x,8) 

Now  substituting  Eqs.  (6.11),  (6.8)  and  (6.6)  into  Eq.  (6.5) 

and  collecting  like  terms  of  order  one  and  <50  yields 
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at  r 


R(x, ©)  (6.12) 


'0  r 


-  R, 


'0  x 


R© 

R2 


*00  Rx 


I  r 


-  Rx?*! 


!e  . 

R2  ^9 


RRX  (  g  ’  ’  <t>0  x  +  iRg  ’  )  COS  © 

Re 

(g*  +  ikg)(cos©  +  sin©) 

at  r  =  R ( x , ©)  (6.13) 


These  are  the  boundary  conditions  for  the  mean  flow  and  the 
cross  flow  for  bodies  of  arbitrary  cross  section.  It  can  be  seen 
that  when  the  body  under  consideration  resumes  its  axisymmetry, 
(i.e.  Rq=0)  Eqs.  (6.13)  and  (6.12)  reduce  to  Eqs.  (2.30)  and 
(2.31)  respectively. 

6.1.2  Pressure  Coefficient 

The  relation  between  the  time  derivative  in  the  wind-fixed 
system,  t’,  and  the  body-fixed  system,  t,  is  provided  by 


0Q 

at  ’ 


00 

at 


g’5(t)r  cos©  +  gi(t)(cos  ©  p  sin  ©  |^) 

(6.14) 


Now  if  Eqs.  (6.10)  and  (6.1)  are  substituted  into  Eq.  (6.14) 
(recall  that  h ( x , t ) = g ( x ) S 0 e 1 k 1 )  and  retaining  only  up  to  linear 
terms  in  <50  yields 
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dfi 

at  ’ 


S  (  t )  <fi .  (  x  ,  r  ,  0 ) 


g  ’  5  ( t )  r  0O  (  x  ,  r  , 9) cos6 


+  i(t)g(^Orcos0 


oesin6) 

/ 


(6.15) 


Substituting  Eqs.  (6.15)  and  (6.11)  into  the  exact  isentrop 
pressure  coefficient  as  given  by  Eq.  (2.32)  and  performing 
binomial  expansion  according  to  the  small  parameter  <S0  ,  we 

Cp  -  Cj  +  <5(t)  C'  (6.16) 

where  C°  =  [S*  -  1]  (6.17) 

oo 

and  Ci  =  -  2S0{(l  +  tf0x)  (*lx-Rg”*Oxcos0)  +  *0  r  (  * ,  r  +g  ’  cos©) 

+  ^o9(^i©-g’sine)  '  ikRg’  ^Oxcos0  + 

+  ikg(0Orcos©  -  sin0^O0)}  at  r  =  R(x,0)  (6.18) 

where  S0  is  now  defined  as 

S 0  =  [1  -  ^2(2^0X  +  *§x  +  *lr  +  p-  *%Q)  }1/Y~1 

The  generalized  forces  will  be  given  now  by 


i  c 


obtain 
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1 


C*  <  '  >  [  (g<  1  >  +  g’<i>  RRx)Rcos©  + 


»2  IT  „  1 


t  j 


1  r  »  f 


j 

j  0 


g( 1 ]R  sinG]  dxdG 
© 


(6. 19) 


Once  the  generalized  forces  are  determined,  the  stability 
derivatives  can  be  obtained  from  Eq.  (2.44),  in  the  similar 
manner  as  to  determine  those  for  axisymmetric  bodies. 


6.2  Mean  Flow  Solution  and  Results 
To  establish  a  robust  unsteady  computational  procedure  for 
asymmetric  bodies  it  is  more  appropriate  to  investigate  first  the 
method  of  solution  of  the  mean  flow  since  there  are  more 
available  data  for  result  verification.  The  mean  flow  problem  is 
governed  by  Eq.  (6.2)  with  the  boundary  condition  Eqs.  (6.4)  and 
(6.12)  . 

If  the  velocities  are  to  be  single-valued  function  the 
dependency  on  the  polar  angle  9  must  occur  through  factors 
cos(n9)  and  sin(n9),  where  n  is  an  integer.  The  elementary 
solutions  in  this  form  that  satisfy  Eq.  (6.2)  were  obtained  by 
Ward  (Ref.  34)  (Chapter  9,  Eq.  (9.3.22))  as 


*(n (x,r,0) 

o 


c os ( n9)  fX  ^ r  1 
sin(n9)J0  ( /?  r  )  n 


\  {(x-<  +  /(x-*)2-/J2r2  )" 


+ 


(x-*-/(  x-<  )2  -/?2  r2  )n  } 


Ldil 

/(x-t)2-/32r2 


df 


(6.20) 
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The  potential  given  by  Eq.  (6.20)  can  be  considered  as  a 
distribution  along  the  body  axis  of  'vortex  multiplets’  of 
strength  f n ( € ) .  For  the  case  n  =  0  Eq.  (6.20)  reduces  to  Eq.  (3.4) 
representing  a  distribution  of  sources  as  presented  in  Chapter  3. 
Since  Eq.  (6.2)  is  linear  a  general  solution  can  be  obtained  by 
superposition  of  elementary  solution,  thus  0o(x,r ,0)  can  be 
expressed 


N 

r  <  n 

0  o  (  X  ,  r  ,  0 )  =  )  ( x  ,  r  ,  e ) 

a  =  0 


(6.21) 


The  functions  f n ( ! )  that  appear  in  Eq.  (6.21)  should  be  obtained 
from  the  boundary  conditions.  But,  before  a  numerical  procedure 
to  determine  them  is  developed,  it  is  convenient  to  estimate  the 
magnitude  of  the  different  terms  in  the  integrand  of  Eq.  (6.20) 
to  see  if  convergency  difficulties  can  arise  during  the  numerical 
solution.  It  can  be  seen  that  the  terms  [  x- f  +  /(  x- (  )  2  ~/3z  r2  ]  n  and 
[ x- C - / ( x- ( ) 2 -p2 r2 ] n  are  of  order  1  since  x  is  defined  between  0 
to  1;  however  the  term  \/{f5r)n  when  evaluated  at  the  body  surface 
is  of  order  l/zn,  where  z  is  the  thickness  of  the  body.  Thus, 
for  thin  bodies  this  term  can  be  very  large  if  the  total  number 
of  elementary  solutions,  N,  is  large.  For  example  for  a  10% 
thick  body  if  N  is  12  those  terms  are  of  order  1012.  Therefore, 
if  the  velocities  have  to  be  of  the  order  of  the  thickness  t 
means  that  the  functions  f n ( f)  must  be  of  order  1  0  ~  1  3  .  This 
implies  that  numerical  difficulties  should  be  expected  if  Eq. 
(6.21)  is  used  to  solve  for  the  mean  flow  potential  <fiQ  ■  To  make 
the  term  1 / z N  small  the  body  thickness  t  has  to  be  large,  which 
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i  s  l  n  contradiction  with  the  assumptions  of  small  disturbances, 


or  N  has  to  be  small.  Since  S  will  depend,  in  general,  on  the 
body  cross  section,  this  latter  condition  can  not  be  guaranteed. 
Thus,  the  application  of  this  method  to  develop  a  robust 
computational  method  should  be  disregarded. 

Thus,  Eq.  (6.2)  should  be  attempted  to  be  solved  by  a  method 
tnat  guarantees  N  to  be  small.  Eq.  (6.21)  can  be  considered, 
with  respect  to  the  dependency  on  0,  as  if  the  mean  flow 
potential  0O  is  fitted  by  cos(nG)  and  sin(nO)  functions  in  the 
interval  0<©<2*r.  Drawing  the  analogy  of  this  method  with  the 
cu r ve- f i t t i n g  method  would  be  like  if  in  the  latter  a  polynomial 
of  h  degree  were  used  to  fit  a  function  s(x).  However  it  is 
known  that  if  a  spline  fit  is  used  the  curve  s(x)  can  be  fitted, 
with  good  accuracy,  by  piecewise  low  order  polynomials.  Thus,  if 
the  concept  of  the  spline  fit  is  used  to  solve  for  Eq.  (6.2)  it 
will  amount  to  divide  the  body  into  interval  along  the  0 
direction,  i.e. 

0  -  0.  <  02  <...<  9J  _  L  <  9,  <  Oj.L  <...<  0N  <  0N  „  !  =  2  rr ,  each 

interval  defined  by  40.  (see  Fig.  36).  On  the  interval 

A0  ,  the  potential  <fr 0  ca  be  expressed,  i.e.  as 

2 

fj./x.r,©)  =  J  90n  j  (x.r.e)  ,  ©j  — ©— ©j  t  ^  ,6.22) 

n  =  0 


where  0‘°,  and  ^0'2  are  given  respectively  by 


x  -  p  r 


f  o .  ,i  LI 


d  f 


(6.23; 


y 1  x- ( ) 2  - p 2  r ‘ 


5  3 


(  1 

>0, 


=  -  cos  0 


'*  /Srl  (x-5)f. 


■xiil  dc 


/(x-«  )2-/82  r2 


(6.24) 


(  2 
*0  , 


=  -  c 


,x-  /Sr 


os20 


(/Sr)2 


(2(x-f)2-/32r2  )  f ,  ,  (5) 
/(x-« ) 2 -02  r2 


<15 


(6.25) 


where  only  cos(n0)  functions  have  been  considered. 

On  each  interval  A0;  a  set  of  M  points  is  distributed  on  the 
body  surface  with  coordinates  (xi  ,ri  ,0^  )  i=l, . . .M,  where 
0j  =  1/2  (0^ +0j  +  j  )  and  ri=R(xi,0J),  such  that  xt<xi  +  1.  Thus,  a 
total  of  M*N  points  on  the  body  surface  are  located  which  are 
called  the  control  points.  From  each  control  point  the 
intersection  of  the  inverse  Mach  cone  with  the  body  axis  is 
determined,  to  obtain  * i + x ,  =  x s -fi r j  with  fx  being  the  body 

apex.  The  interval  to  J  is  called  a  panel.  On  each 

panel  the  strengths  of  the  sources  (f0>J),  dipoles  (f^ j)  and 
quadrupoles  ( f2  )  is  assumed  to  vary  linearly,  i.e. 


fo, j(«) 

=  a> 

+  b  . 

J 

fl.  j(») 

=  c. 

.  j* 

-r  d  ! 

1 ' J 

5<5<5l 

tl  (6.26) 

f2,  ,  (<) 

=  e  i 

. 

1 

+  J 

If  the  condition 

that 

the 

functions  f0iJ»  fx  t  j  , 

and  f  2 (  j  must 

continuous  between  each  two  panels  is  imposed,  the  constants  btj 
can  be  expressed  as  a  function  of  a t  ^  ,  b x  ,  and  d x  }  of  c t  ^  , 

d,  andi  ofe  ,  ■  The  constants  b,  ,  ,  d .  ,  and  i. 

can  be  shown  to  be  zero  when  the  boundary  condition  at  the  Mach 
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wave  emanating  from  the  apex  is  satisfied.  Thus  the  potential  at 
the  control  point  xn ,  rn ,  Q,  can  be  shown  to  be  after  integration 


n  _ 

\  ~  {a.  ,  j  )cosh‘  1 


1 


-  /(x„-S  )2~P2rl  + 


J  L 


/Jr  cosh'  :  ^ — — 

Prn. 


cos 0^  +  ej | j 


Pr_a 

8 


'lf;(xn-t)^(Xn-(^^  * 

-3*-  ^rn 


-  (' 


11^. 


xn-«  )  +  /(x„-€  )2-/J2  r2 


)3} 


_££a_ 


(xn  -«)  +  /(  xn-«)2'/J2r2 


-  (■ 


(xn-i)+/(xn-f )2-^2r2 


)  j  cos ( 20J ) j 


i  ♦  1  .  J 

i  .  J 

(6.27) 


The  velocities  in  the  x,r,0  directions  can  be  obtained  by 
applying  the  derivatives  with  respect  to  x,r,  3nd  0  respectively 
Thus,  the  mean  flow  potential  and  velocities  can  be  expressed  as 
function  of  3M*N  unknown  constants  al  ,  cl:J  and  ei,j  i  =  l»---M, 
,j  =  l,...N.  If  the  conditions  that  the  potential  as  well  as  the 

tangential  velocity  in  the  0  direction  must  remain  continuous 
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between  adjacent  intervals  i.e.  A©j  to  A© ; +  ^  are  imposed,  a  set 
of  2M*N  equations  relating  a(i  j  ,  ci  ,  j  and  ei  j  are  obtained. 

These  equations  can  be  used  to  express,  i.e.  the  c;  ’  s  and  the 
eiJ’s  as  functions  of  the  ai _  ’ s .  Now,  by  satisfying  the 
tangency  condition  at  the  control  points  the  a,  j’s  can  be 
determined.  Once  the  a;  ^  ’ s  are  known,  the  velocities  at  the 
control  points  can  be  evaluated  and  from  Eq.  (6.17)  the  pressure 
coefficient. 

It  should  be  noted  that  the  number  of  elementary  solutions 
chosen,  three  in  the  above  description,  can  be  altered.  For 
example,  if  instead  of  imposing  the  continuity  in  the  tangential 
velocity  in  the  ©  direction,  only  the  potential  is  required  to  be 
continuous  between  adjacent  intervals,  only  two  elementary 
solutions,  i.e.  sources  and  dipoles  need  to  be  considered  on  each 
interval.  Also  any  other  set  of  vortex  multiplets,  instead  of 
sources,  doublets  and  quadrupoles,  can  be  used  and  the  type  of 
vortex  multiplets  used  can  change  from  one  interval,  A Qjt  to  the 
next  AQj + ! .  This  method  can  be  called  spline  panel  method. 

The  simplest  application  of  the  method  of  solution  described 
above  is  if  only  one  elementary  solution  on  each  interval  is 
considered.  In  that  case,  the  potential  a r. d  its  derivatives  in 
the  ©  direction  are  not  continuous  between  adjacent  intervals. 

The  best  choice  of  elementary  solutions  is  the  source 
distribution,  since  it  is  known  to  be  the  exact  solution  of  the 
mean  flow  potential  for  bodies  of  revolution.  Thus,  the  use  of 
source  singularities  of  different  strength  on  each  interval  A©, 
amounts  to  treating  locally  the  body  as  a  body  of  revolution  with 
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radius  R(x,©  ).  Therefore,  some  restrictions  regarding  the 
arbitrariness  of  the  body  cross  section  at  which  this  first 
approximation  can  be  applied,  should  be  expected. 

In  the  Fig.  37  the  mean  flow  pressure  for  a  conical  body  at 
Mach  number  M  =2.0  whose  cross  section  varies  as  cos©  is 

oo 

compared,  at  x=l,  with  Devan’s  results  (Ref.  35),  and  USSAERO 
results.  It  can  be  seen  that  except  in  the  region  8* <©<30* ,  the 
agreement  seems  to  be  very  good  with  both  methods.  The  deviation 
in  that  region  we  believe  is  due  to  the  use  of  only  one 
elementary  solution  on  each  interval.  However  it  is  seen  that 
this  small  deviation  does  not  affect  in  the  prediction  of  the 
axial  and  normal  forces  when  compared  to  the  other  two  methods. 
The  small  difference  between  Devan’s  results  and  USSAERO  results 
can  be  due  to  the  different  method  used  to  solve  the  differential 
equation.  The  first  used  the  finite  difference  method  while  the 
second  used  the  panel  surface  method.  The  computer  time  taken  by 
these  two  methods  is  at  least  one  order  of  magnitude  larger  than 
the  present  method. 

In  the  Fig.  38  the  mean  flow  pressure  for  a  conical  body  of 
circular  cross  section  for  -n/2ZQZn/2  and  elliptical,  with 
ellipticity  ratio  1/2,  for  rr/2<Q<3/2rr  at  Mach  number  Mw>=2.0  is 
compared  with  Devan’s  (Ref.  35)  and  USSAERO  results.  It  can  be 
seen  that  a  slight  deviation  occurs  overall  the  values  of  ©  when 
compared  to  the  other  two  methods  and  the  normal  force 
coefficient  has  a  difference  of  about  6%.  However  the  axial 
force  coefficient  is  very  well  predicted.  The  reason  for  the 
discrepancy  is  that  an  ellipticity  ratio  of  1/2  is,  too  low  for 
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the  present  method  to  produce  accurate  results.  As  will  be  shown 
later  for  ellipticity  ratio  of  0.75  to  1.3,  the  present  method 
gives  good  agreement  when  compared  to  USSAERO  results. 


6.3  Unsteady  Flow  Solution  and  Results 
The  unsteady  flow  problem  is  governed  by  Eq.  (6.3)  with  the 
boundary  conditions  Eqs.  (6.4)  and  (6.13).  Like  for  bodies  of 
revolution.  Chapter  4,  only  the  homogeneous  unsteady  wave 
equation  will  be  solved.  The  term  G3(g,0o)  can  be  taken  into 
account  in  a  similar  manner  to  that  described  for  bodies  of 
revolution  Eq.  (4.1). 

To  extend  the  method  described  in  the  previous  section  to 
unsteady  flow  it  is  necessary  first  to  obtain  a  general 
elementary  solution,  equivalent  to  Eq.  (6.20)  for  the  mean  flow, 
for  the  unsteady  flow.  This  solution  can  be  assumed  in  the  form 


(  n 

(x,r,9): 


cos(ne) fX  .  -ip(x-l) 

s  in  (  n9)  J  fn'f)e 


cos(  <V(x-f)2-,g2  r2  ) 
/( x-< ) 2 ~p2 r2 


m(cosh'  1  — — )  ;  n)  dt 


(6.28) 


where 


e~i m(x-«) 


cos(  it /( x  -  l )  2  -  p2  r2  ) 
/(x-t)2-/?2  r2 


represents  the  kernel 


for  a  distribution  of  unsteady  sources  along  the  body  axis,  and 
it  is  known  to  satisfy  the  equation  (see  Garrick  (Ref.  26)) 


(  1  “M2  )  +  d> ,  +  —  a J 

v  ooy  'lxx  ^1  rr  r  "lr 

-  ^-=-  2ikM2  <j> ,  +  k2M2  <j> .  -  Q  when  n  =  0.  (6.29) 

J"  fc  »  OO  X  X  OO  x 


Thus,  in  Eq.  (6.28)  m  represents  a  function  of  x,r,€  and  n 
to  modify  the  source  kernel  so  it  satisfies  Eq.  (6.29)  for  n 
different  of  zero.  Then,  m  must  satisfy  m(x,?,r;  0)=1.  To 

determine  m  we  first  make  the  following  transformation  in  Eq. 
(6.28) 


f  =  x  -  pr  cosh  a 

di  =  -  pr  sinh  a  d  a 

at  i  =  0 

a  =  cosh'  1 

pr 

and  at  «  =  x-fiv 

a  =  0 

Eq.  (6.28)  then  becomes 


<  n 

4>x  (x,  r,0) 


cos ( n©)  r 
s in ( n©) j 


cosh 


>(£-) 


Pr 


F  (  x-/Sr  cosha )  e 


-i/i|Srcosh<r 


cos ( f^rs inha)  m ( a ; n )  da 


(6.30) 


To  determine  m(a;n)  the  derivatives  with  respect  to  x  and  r 
in  Eq.  (6.30)  are  determined  and  substituted  into  Eq.  (6.29).  In 
this  way  it  can  be  shown  that  Eq.  (6.30)  satisfies  Eq.  (6.29)  if 
m’’-n2m=0.  Therefore,  the  solution  to  this  equation  together 
with  the  condition  m(a;0)=l,  implies  that  m  should  be 
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m  ( (7  ;  n  )  =  cos  h  (  na )  .  Thus,  the  generalized  elementary  solution  of 
Eq.  (6.29)  is  given  by 


(  n 

<t>x  (  X  ,  r  ,  0) 


cos ( n©) 
s i n ( n©) 


•cosh 


i  (*_) 


,  ,  .  -i^/Srcoshcr 

(  x  -  p  r  c  o  s  h  <7 )  e 


cos(  r'/Srsinhff)  cosh(nff)  da 


(6.31) 


Now,  the  extension  of  the  spline  panel  method  to  unsteady 
flow  can  be  done  by  letting 

2 

4 i , j ( x , r , ©)  =  ^  *i,j(x»r,©)  ,  ©j  <  ©  <  ©j + x  (6.32) 

n  =  0 


where  <fi  l  _  ■  means  the  unsteady  flow  potential  at  the  interval  j. 

(x,r ,©)  represents  a  distribution  of  unsteady  sources  along 
the  body  axis,  1  of  unsteady  doublets  (it  can  be  shown  that 
is  equal  to  that  given  by  the  first  term  in  Eq.  (4.1)  for 
bodies  of  revolution  by  integrating  by  parts),  and  <j>\z  of 
unsteady  quadrupoles.  The  Harmonic  Gradient  concept  needs  now  to 
be  applied  to  the  strengths  F0  j  .  ( I )  ,  and  F2  (  j  (  *  )  in  a 

similar  manner  as  was  done  in  Chapter  4,  Eq.  (4.6)  for  FX(S). 

The  integration  of  0 ,  and  2  will  not  be  carried  out 
here;  however  it  is  similar  to  that  shown  in  Appendix  B  for 
bodies  of  revolution  (which  is  the  same  as  for  <j>\  1  )  .  Once  the 
integration  is  done  the  method  of  solution  to  determine  ,  and 
the  velocities  x  ,  r  and  l/r#1Q,  at  the  control  points  is  the 
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same  as  the  one  described  in  the  previous  section  for  the  mean 


flow. 

Since  for  the  mean  flow  numerical  results  were  determined 
using  only  one  elementary  solution  on  each  interval  AQj  , 
consistently  the  same  will  be  done  for  the  unsteady  flow.  Thus, 
on  each  interval  a  distribution  of  unsteady  doublets  of  strength 
Fx  (5)  is  distributed.  Like  for  bodies  of  revolution  the 
doublet  strength  is  discretized  along  the  body  axis  and  the 
Harmonic  Gradient  concept  is  applied. 

For  the  unsteady  flow  a  set  of  conical  bodies  of  elliptic 
cross  section  and  ellipticity  ratio  a/b  from  0.75  to  1.3  have 
been  investigated,  at  Mach  number  M^-3. 0,  and  compared,  whenever 
possible  with  USSAERO  results. 

In  the  Figs.  39  and  40  the  static  normal  force  and  moment 
coefficients  are  presented  versus  a/b  and  compared  to  those 
obtained  from  the  USSAERO  code.  It  can  be  seen  that  good 
agreement  is  obtained  for  both  coefficients  within  the  range 
studied. 

In  the  Fig.  41  the  dynamic  normal  force  and  moment 
coefficients  are  presented  versus  a/b.  No  comparison  with  other 
theories  is  made  since,  apparently,  this  is  the  first  method  able 
to  compute  the  unsteady  aerodynamics,  for  bodies  of  arbitrary 
cross  section,  in  supersonic  flow.  At  the  present  time  no 
experimental  data  for  the  unsteady  coefficients  could  be  found. 
However  since  the  static  coefficients,  within  the  range  studied, 
are  in  good  agreement  with  the  USSAERO  results,  it  is  expected 
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that  the  unsteady  coefficients  should  also  be  well  predicted  by 
the  present  method. 

In  the  Figs.  42,  43  and  44  the  real  and  imaginary  part  of 
the  pressure  coefficient  for  the  elliptic  cones  in  rigid,  first 
bending,  and  second  bending  mode  respectively  are  presented  at 
0=0,  and  reduced  frequency  k=1.0,  along  the  x  axis.  In  general 
we  can  see  that  the  imaginary  part  of  the  pressure  coefficient 
slightly  changes  with  the  ellipticity  ratio,  and  that  when  a/b 
approaches  one,  the  present  results  approach,  as  expected,  those 
obtained  for  the  body  of  revolution. 

In  the  Figs.  45,  46  and  47,  the  real  and  imaginary  part  of 
the  pressure  coefficient,  for  rigid,  first  bending  and  second 
bending  mode  respectively  are  presented  versus  ©  at  x=l.  The 
real  part  of  the  pressure  coefficient  for  the  case  a/b=0.75 
deserves  some  attention.  To  investigate  this  behavior  the 
pressure  coefficient  for  an  elliptic  cone  (a/b=.75)  and  a 
circular  cone  at  steady  angle  of  attack  are  compared  to  the 
pressure  coefficient  for  a  flat  plate,  obtained  from  the  slender 
wing  theory  (Ref.  5)  along  the  y  axis  in  Fig.  48.  Although  for 
the  present  time  lower  values  of  a/b  have  not  been  obtained,  it 
can  be  seen  that  the  trend  in  pressure  of  the  present  method  is 
the  same  as  that  obtained  from  the  slender  wing  theory  except  in 
the  region  where  y/b  approaches  one,  where,  for  a  flat  plate,  the 
pressure  becomes  singular. 

It  is  interesting  to  notice  that  when  a/b  decreases  the  body 
becomes,  as  expected,  more  stable  statically  and  dynamically. 
However,  from  the  structural  point  of  view  we  would  say  that  the 
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body  becomes  more  flexible  and  it  might  turn  out  to  be  more 
likely  to  flutter.  Since  the  present  approach  has  no  limitations 
regarding  frequency  or  mode  shape  it  can  be  an  excellent  tool  to 
estimate  the  advantages,  from  the  stability  point  of  view,  and 
the  disadvantages,  from  the  aeroelastic  point  of  view,  of  letting 
the  ellipticity  ratio  to  decrease  and,  thus,  to  obtain  an 
optimized  cross  section  body  shape. 
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CHAPTER  7 


BUNDLED  TRIPLET  METHOD 

7.1  New  Development  of  BTM 

As  mentioned  in  Chapter  6  using  Fourier  representation  of  an 
arbitrary  body  shape  would  introduce  a  formulation  involving 
higher-order  kernel  function,  which  makes  the  determinant 
evaluation  very  difficult.  On  the  other  hand,  the  Spline d- Panel 
Method  used  in  the  same  chapter  has  shown  limitations  on  the 
degree  of  asymmetry  for  bodies,  which  indicates  that  this  method 
may  not  be  general  enough  to  handle  arbitrary  body-wing 
combinations.  The  proper  treatment  of  the  body-wing  interference 
mandates  to  generalize  the  Spline-Panel  method  to  the  Bundle 
Triplet  Method  (BTM)  in  which  at  least  two  substantial  features 
have  been  added.  First,  to  keep  the  kernel  function  in  low 
order,  multiple  lines  of  low  order  singularities,  namely  line 
sources  and  line  doublets,  are  employed.  Second,  by  using  a 
least  square  procedure,  the  BTM  can  account  for  the 
c i rcumf e r enc i a  1  influence  between  panel  elements,  thus  the  method 
is  sufficient  to  treat  bodies  of  arbitrary  cross-sections.  As 
can  be  seen  in  the  later  computed  cases,  the  results  obtained  by 
the  BTM  are  all  in  good  agreement  with  existing  theories  and 
measured  data. 

7.2  Formulation 

As  shown  in  Figs.  49  and  50  the  body  cross  section  is 
divided  into  'M’  intervals.  Each  interval  contains  a  finite 
sector  40  ,  where  40  =  0„  ^  .  -  0_ ,  m  =  1 ,  ...  M,  and  0,  =  0, 

m  1  n  in+i  m  *  1  *  1  • 
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=  2/r.  A  three-dimensional  "pie”  shape  is  defined  by  a 
portion  of  the  body  within  the  sector  A3  m  whose  vertex  lies  along 
the  body  axis.  A  line  distribution  of  sources  and  another  of 
doublets  are  superposed  along  the  body  axis  for  each  pie  segment. 
Within  each  pie  segment  an  integral  solution  can  be  expressed  as 
based  on  the  superposition  principle, 


x,-£ri 

*1 ^X1  - r I  - )  =  J  Fa(f)B(xl  -  t,  fir i  ) d  € 
0 


+  J  G0(«)^r-B(xt  -  f,  ^r,  )  df  •  cos*,  (7.1) 


<3r 


+  J  Hm  (  -  i,  Arjdl-sin  »t 

0 


dr 


where  f  x  t  ,  r4,  8 m)  represents  a  typical  control  point  lying  on 


the  body  surface  r4 
kernel  function  B(xi 

B  (  x,  -  « ,  firi  ) 


=  R(x,  )  and  8m  =  ( Sm  +  8  B  +  x  )/2. 
-  t,  /Srl  )  can  be  expressed  as 

_  -i^(x,-«)  cos ( xRn ) 

R0 


The 


where  RQ  =  {(xi  -  *)2  -  pz  T ^ / z 


(7.2) 


M  =  kM*/^  ,  X  -  kM J  ^ 


Fm ( f )  denotes  the  source  strength  distribution  and  Gm ( f  )  and 
Ha ( f )  the  doublet  strengths  distributions  in  the  ml h  interval. 


The  fact  that  the  doublet  distribution  has  two  unknown 


functions  can  be  interpreted  as  that  the  doublet  strength  as  well 
as  the  direction  of  the  doublet  axis  are  required  to  be 
determined.  It  is  clear  now,  our  triplet  model  is  a  "generalized 
triplet"  in  that  the  singularity  strength  functions  are  different 
from  each  other.  In  the  case  of  Fm  =  Gm  -  Hm ,  the  generalized 
triplet  model  then  reduces  to  the  regular  triplet  model  (see  Ref. 
36)  . 


Furthermore,  along  the  body  axis,  each  of  the  triplet  lines 
is  discretized  by  N  segments  where  AS ^  _ m  =  S  + x , m  -  f  m  for  a 
given  pie  sector  A 8^  (See  Figs.  49  and  50).  Also  along  the  body 
axis,  different  source  and  doublet  strength  functions  are 
discretized  as  Fm  -*  (  t )  ,  Gm  j  (  S  )  and  Hm  J  (  S )  within  AS-  m  .  The 
velocity  potential  within  this  sector  can  then  be  expressed  as 

*  j  *  i 

*i  (x,  ,rt  .O  =  J  J  FJ|B(f)  B(xt  -  5,  Arjdf 


j  =  l 
1 

+  ^  I  GJ  ,  rn  (  1 )  §7  B  (xi  ~  Pri  )d«  *cos^  (7.3) 

J  =  1 


i 


j  ♦  i 

H  , 


(  ? 


S ,  fir  x  ) ds • sinff m 


7.3  Harmonic  Gradient  Model 

To  achieve  computation  accuracy  and  effectiveness  in  the 
high-reduced  frequency  and/or  the  low  Mach  number  domains,  it  is 
essential  to  render  the  doublet  solution  uniformly  valid  in  the 


complete  k  and  M  domains.  This  can  be  simply  achieved  by 
modeling  the  doublet  solution  to  maintain  spatially  harmonic  in 
the  mean  flow  direction.  In  so  doing,  the  panel  size  A f J  is 
regulated  and  is  made  compatible  to  the  wave  number  generated  by 
the  body  oscillation.  This  is  known  as  the  Harmonic  Gradient 
(H-G)  model  introduced  in  Ref.  4. 

Here,  the  H-G  model  employed  for  wings  (Ref.  4)  is  adopted 
to  combine  with  the  Bundle  Triplet  solutions.  Application  of  the 
H-G  model  to  Eq.  (7.3)  amounts  to  stating  that  in  the  interval 

5  ;  .  m  V  ^  j  +  1  ,  m  1 

-i/i? 


=  aj , m ( 1  - 


for  the  source  solution,  and 
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=  (b4,.«  +  d4l.)' 


•i/»(xi  -f ) 
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—  =  vC 


€  +  h 


j  ,  m 


)e-i/i(Xi-0 


(7.4b) 


(7.4c) 


for  the  doublet  solutions. 

Note  that  Eqs.  (7.4b)  and  (7.4c)  are  the  direct  application 
of  the  H-G  model  whereas  Eq.  (7.4a)  is  a  new  extension  of  the  H-G 
model  to  the  source  solution.  It  can  be  shown  that  both 
constants  d .  m  and  hj  n  can  be  expressed  in  terms  of  bj  a  and 
cj  a  respectively  when  a  continuity  condition  in  Gn  and  Hm  is 
imposed  between  the  adjacent  segments.  Introducing  Eqs. 

' 7 . 4  a , b , c !  to,  and  further  discretization  of  the  kernel  integrals 
of  Eq.  '7.3;  result  in 
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3  +  1 


(7.5) 
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j  =  l 


+  (bJ)aCOS(?t 


i  s in^.  )  D j  ') 


j  +  l  j  +  i 

where  S  j  and  D;  represent  the  induced  potentials  per  unit 
source  and  doublet  strengths  respectively  by  the  jir>  segment  of 

the  triplet  line  at  the  point  (xt  ,r.  ,0  );  they  are  given  by 


j  +  i 


j  ♦  i 


-  e 


iM 


- 1  co"XKo  d* 


and 


j  +  l 


'  j  *  i 


dr, 


S  (x;  -  £  ,  /Sr  l  )  d£ 


where 


x  ,-£ 


S(x,  -  I,  fr,)  -  f  ZZLi d, 

^r. 


(7.6) 


(7.7) 


Evaluations  of  DjJ+1  as  well  as  its  derivatives,  with 
respect  to  xt  and  rt  ,  are  identical  to  those  shown  for  bodies  of 

j  ♦  l 

revolution  in  Appendix  B.  For  detailed  evaluations  of  Sj  and 
its  derivatives,  with  respect  to  and  rt  ,  one  is  referred  to 

Appendix  D.  There  it  is  shown  that  after  a  transformation  these 

integrals  can  easily  be  evaluated  numerically  by  Gaussian 

quadrature . 

The  domain  of  influence  is  defined  by  a  pie  segment  bounded 
by  the  sector  and  the  inversed  Mach  cone  intersection  with 
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respect  to  P , .  Hence  the  effective  triplet  line  is  confined  to 
the  part  from  the  body  apex  to  the  iLh  segment  where  the 

i  +  1  j  +  1 

discretized  kernel  integrals  Sj  and  D.  represent  the 
induced  potential  per  unit  source  and  doublet  strengths  by  the 


i  t  h 


segment  of  the  triplet  line.  Moreover,  it  should  be  noted 


that  the  steady  mean  flow  potential  </>0  at  the  point  (x,  ,  rL  ,  %  ) 

j  ♦  1  J  ♦  1 

can  be  obtained  from  Eq.  (7.5)  by  letting  k=0  in  Sj  and  Dj 

7.4  Least  Square  Procedure 

For  convenience,  a  surface  panel  is  defined  as  a  body 
surface  element  whose  area  is  bounded  by  R  •  and  hfj  .  A 

control  point  is  located  at  (xf  ,Ri  (x4  ,  9m)  ,  9m)  ,  the  central  point 
of  each  panel,  where  i  =  l,...N,  m=l,...M,  and  Sm = ( 9m  +  9m  +  1 )  / 2 . 

To  obtain  the  potential  values  of  <f>l  _  m  at  the  control  point,  one 
of  the  effective  ways  is  to  relate  the  unknowns  a^  f  m  ,  bj _ m  and 
Cj  m  to  the  cross  flow  potential  by  means  of  the  method  of 

weighted  least  squares.  (For  simplicity,  we  have  dropped  the 
subscript  "1"  in  hereafter.) 

To  apply  the  least  square  method  one  considers  the  induced 
potential  by  the  singularities  distributions  at  the  h  pie 
segment  at  the  points  ( x ,  , r ,  ,  3 m _ l  ) *  and  ( x ,  , r t  , 9m „ l  )  ,  i.e. 

*1  '  r.  '  *m-  1  >  =  Z  aj.mSj  ^  <bj  ,  m  cos^m-  l  +Cj  )D' 

j  =  1 


and 


(7.7) 


'l  U,  ,r,  )  =  2  a  J  ,  m  ^  j  l+  (bj  ,  o»COS*m.  1+CJ  ,  „  S  i  n  <?  m  ,  1  )  D  J 

)  =  1 


It  should  be  noted  that  for  this  point  rl  =R(xt  , 9m _  x  )  and  for  the 
next  r,  =  R(x,  ,  j  )  . 
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j 


Let  h  be  the  exact  value  of  the  potential  at  the  body 
surface,  the  source-  and  doub 1 e t -s t r en g t h  solutions,  due  to  the 
lease  square  approach,  arc  the  ones  that  minimize  the  residual 


I  =  k 


m  +  1 
X 

'h  =  m  - 


Wh{[  2 


j  +  l 

,S,  +( b 


j  =  i 


mCOS*h+Cj 


asin*h 


]-( 


}2  (7.8) 


J  *  1  J  1 

where  Sj  and  Dj  are  now  evaluated  at  the  control  point 
( x i , R t ( x , , & h ) 9 h ) ;  ^  h  represents  the  exact  value  of  the 

potential  at  the  ith  point  and  Wh  are  the  introduced  weight 

functions.  Note  that  Eq.  (7.8)  is  a  least  square  procedure  set 

up  for  three  adjacent  panels,  namely  from  ( m- 1 )  to  (m+1).  Thus, 

imposing 


61 _  61 

<3a,  <3b  j 

i  .  m  i  .  n 


(7.9) 


yield , 


"z1  Wh  {X  3j  ,  mS^  +  1  +(b0  ,  mcos#h  +Cj  ,  nsin^  )Dj  l-*llh}s|  '  =0 

h  =  m  -  L  j  =  1 


“x1  Wh{Xaj>msJ  1  +(bj  ,  Bcos#h  +Cj  _  msin*h  )  Dj  1 ilh}D‘  ‘cos^  =0 

h  =  m  -  1  j  =  1 


■x‘  Wh(2aj|fflsJ  1  +  (  bo  ,  Dcos«h  +Cj  _  nsin*h  )  dJ 

h  =  m  -  l  j  =  1 


.  h  )  D!  =0 


(7.10) 


These  equations  can  be  expressed  in  matrix  form  as 


[A] 


j  ,  B  1  1-1 

raj  ,  o  i 

ril  i  ,  m  -  1 

J..}  +  Z  CY]- 

[H], 

{*L  .1 

j  j  -  I 

j  ,  m 

cj  ,  m 

^  i  ,  m  ♦  1 

i  A  ]  i  i  [  7  1  j  ,  i 

and  [ H ]  t 

are 

given  in 

(7.11) 
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We  can  now  relate  the  unknowns  a  m  , 

i  ,  m  7 


where  [U],  [V]  and  [W]  are  matrices  of  dimension  ( NxM ) x ( 3xNxM ) 

containing  the  velocity  influence  coefficients  in  the  x,  r  and  9 

direction.  Now,  the  velocities  ^ —  and  - —  —  in  the 

dx,  drl  rx  d9m 

tnagency  condition,  Eq.  (6.13),  can  be  replaced  by  the  cross  flow 
potential  <f>  i  _  m  through  Eqs.  (7.13-15).  Consequently,  m }  is 

the  unknown  to  be  solved  in  Eq.  (6.13) 

In  passing,  we  note  that  the  evaluation  of  the  steady  mean 
flow  potential  follows  the  same  procedure  as  described  above, 


7  1 


J  +  1  J  +  1 

except  that  the  kernel  integrals  Sj  and  D  are  to  be 
replaced  by  their  steady  counterparts. 

7.5  Panel  Flutter 

To  verify  the  solutions  obtained  in  the  high  reduced- 
frequency  range,  one  is  required  to  apply  the  present  BTM  to  the 
supersonic  panel  flutter  problem  for  cylindrical  shells. 
Previously,  Dowell,  and  Widnall  (Ref.  37)  used  Laplace  transfer 
technique  to  study  this  problem,  whereas  Platzer  et  al .  (Ref.  38) 
used  the  Linearized  Method  of  Characteristics  ( LMOC ) .  The 
objective  of  these  studies  is  to  evaluate  the  generalized  forces 
acting  on  the  vibrating  cylindrical  panel  of  an  otherwise  rigid, 
infinitely  long,  cylindrical  shell. 

Let  L  be  the  length  of  the  cylindrical  panel  i.e. 

(0  <  x  <  L).  Along  the  cylindrical  shell,  the  oncoming  flow 
upstream  of  the  panel  is  uniformly  supersonic.  The  cylindrical 
panel  is  allowed  to  perform  smal 1 -amp  1 i tude  harmonic  oscillations 
(see  Fig .  55 ) . 

For  inviscid  analysis  of  panel  flutter,  Eq.  (6.3)  (G.3  =  0) 

is  the  commonly  used  governing  equation.  Or.  the  mean  surface  of 
the  cylindrical  panel,  r=R,  the  tangency  condition  reads: 

4>x  r  =  *kh  +  hx  ,  0  <  x  <  L  (7.16) 

=  0  ,  x  <  0 

where  h  is  the  mode  shape  of  the  vibrating  panel,  defined  as 

h(x,0)  =  s  i  n  (  j  nx. )  cos  (  n  9  )  (7.17) 
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The  unsteady  pressure  coefficient  Cp  on  the  panel  surface  is 
simply 

C  p  =  -2 [ 0  x  x  +  ikrfj  (7.18) 

and  the  generalized  aerodynamic  forces  Q  reads 
2  n  L 

'  1  (  i  ) 

Cp  s  i  n  (  j  rrx  )  cos  (  n8  )  R  d0  dx  (7.19) 

3 

i(i> 

where  Cp  is  the  unsteady  aerodynamic  pressure  due  to  the  mode 

sin(iTx)cos(n»)  and  C  =  2n  for  n  =  0  and  tt  for  n  *  0. 

Bundled  triplets  are  placed  along  the  x-axis  in  the  interval 

-/SR  <  i  <  L-/SR.  The  panel  elements  and  control  points  between 

0  <  x  <  L  and  0  <  9  <  2n  on  the  mean  surface  are  distributed 

according  to  the  cosine  law  in  both  x  and  9  directions. 

Combining  Eqs.  (7.14),  (7.16)  and  (7.17)  results  in  the  following 

tangency  condition  evaluated  at  the  control  points. 

(  V]  [  LS  ]  [<fii  '  m  }  =  {  (  ik  s  in(  jirxi  )  +  jn  cos  (  j  /rxi  )  )  cos  (  n^  )  }  (7.20) 

The  potential  values  4>l:m  can  now  be  obtained  from  the  above 
equation.  Thus,  the  unsteady  pressure  and  the  generalized  forces 
can  be  evaluated  according  to  Eqs.  (7.18)  and  (7.19). 

In  the  present  formulation,  the  circumferential  mode  shape 
(i.e.,  the  cos(nS)  factor)  is  retained  throughout  the  analysis, 
whereas  Refs.  37  and  38,  cos(nS)  is  factorized  out  in  their 
formulation.  As  a  result,  the  number  of  control  points  required 
in  the  circumferential  direction  increases  rapidly  with  an 
increase  in  n,  although  the  computation  time  is  still 


1 

2C 
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comparatively  rapid.  It  is  essential  to  point  out  that  the 
present  method  is  a  full  three  dimensional  one  which  can  treat  a 
much  wider  scope  of  problems  than  those  in  Refs.  37  and  38. 
Feasible  applications  in  panel  flutter  using  the  present  method 
include  closed  (pointed)  bodies  of  circular  or  noncircular  cross 
sections  and  vibrating  panels  of  non-harmonic  mode  shapes  (i.e., 
the  case  of  partially  oscillating  shells). 

7.6  Results  and  Discussion 

To  validate  the  present  method,  numerical  examples  in  terms 
of  steady  pressures,  stability  derivatives,  and  generalized 
forces  are  presented  for  asymmetric  bodies  and  for  cylindrical 
panels . 

7.6.1  Asymmetric  Bodies 

The  steady  pressure  distributions  for  three  different 
asymmetric  conical  bodies  are  presented  in  Figs.  51,  52  and  53  at 
the  same  freestream  Mach  number  M  =  2.0.  Since  the  flowfields 

o© 

are  conical,  only  the  circumferential  Cps  are  presented  for  these 
cases.  Asymmetric  configurations  as  shown  in  Figs.  51  and  52  are 
placed  at  mean  angle  of  attack  a  -  0 ,  whereas  the  elliptic  cone 
in  Figure  4  is  at  a  =  5°.  Since  the  flow  is  symmetrical  with 
respect  to  the  meridian  plane,  pressures  on  half  of  the  body  are 
presented  (0  <  9  <  180°).  Because  of  the  steeper  variation  in 
the  given  body  curvature  (due  to  cos  30),  a  bundle  of  36  triplet 
lines  is  distributed  in  equal  circumferential  intervals  for  a 
full  body  in  Fig.  51.  The  geometries  in  Figs.  52  and  53  are  less 
complicated;  only  18  triplet  lines  were  used  for  full  bodies  in 
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both  cases  . 


It  can  be  seen  that  the  present  computed  results  are 


in  good  agreement  with  Devan’s  results  using  the  finite 
difference  method  (Ref.  35)  and  the  USSAERO  results  (Ref.  2). 
Typically,  we  use  40  panels  in  the  circumferential  direction  in 
the  USSAERO  program  for  obtaining  results  in  Figs.  51,  52  and  53. 

In  the  present  method,  5  segments  are  used  in  the  x 
direction.  This  amounts  to  a  total  of  100~200  panels  and  control 
points  to  be  evaluated.  Because  the  evaluation  scheme  of  the 
present  kernel  is  much  simpler,  the  CPU  time  required  is  about 
one  tenth  of  that  needed  in  the  USSAERO  code. 

Fig.  54  presents  the  static  (Figs.  54a,  54b)  and  dynamic 
(Figs.  54c,  54d)  normal  force  and  moment  coefficients  at  various 
Mach  numbers  for  a  family  of  elliptic  cones  placed  at  zero  mean 
angle  of  attack.  The  ellipticity  ratio  as  defined  by  a/b ,  ranges 
from  .75  to  1.3.  Figs.  54a  and  54b  compare  the  present  results 
with  those  computed  by  USSAERO.  There  appears  to  be  an  increase 
in  discrepancies  as  the  elliptic  cone  becomes  more  wing-like, 
i.e. ,  a/b  >  1.  When  the  ratio  approaches  one,  all  results 
converge  to  the  results  for  a  circular  cone  as  obtained  in 
Chapter  4 . 

Similarly,  the  computed  results  for  the  out-of-phase  normal 
forces  and  moment  coefficients  (Figs.  54c  and  54d)  also  converge 
to  those  obtained  in  Chapter  4  when  a/b  approaches  one  as 
expected.  In  passing,  we  note  that  little  Mach  dependency  is 
observed  for  the  static  moments  (about  xG  =  0),  whereas  the 
damping  moment  decreases  rapidly  with  increasing  Mach  number. 
Furthermore,  all  results  confirm  the  expected  trend  that  both  the 
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static  and  the  damping  moments  increase  with  decreasing 
ellipticity  ratio,  i.e.,  as  the  body  becomes  more  wing-like. 

7.6.2  _ Cylindrical  Panel  Flutter 

In  order  to  validate  the  present  method  in  the  high- 
frequency  domain,  supersonic  cylindrical  flutter  cases  (Fig.  55) 
are  selected  for  comparison  with  existing  theories  (Refs  37  and 
38)  . 

A  bundle  of  triplet  lines  are  arranged  according  to  a  cosine 
distribution  both  in  the  circumferential  and  in  the  axial 
directions  (0  <  8  <  2rr  and  0  <  x  <  1).  The  cylindrical  panel  is 
first  evenly  divided  into  n  intervals  in  the  circumferential 
direction,  say  n=5.  Within  each  interval,  eight  control  points 
are  used,  given  a  total  of  40  control  points.  In  the  axial 
direction,  25  points  are  used  for  all  n’s. 

The  real  and  the  imaginary  parts  of  the  generalized  forces 
on  the  cylindrical  panel  are  presented  in  Figs.  56  and  57  for  the 
freestream  Mach  number  =  /2  and  the  reduced  frequencies  k=0 
and  1.0.  The  generalized  forces  are  computed  based  on  Eq. 

(7.19). 

It  is  seen  that  good  agreements  are  found  between  the 
present  results  and  those  of  Dowell  and  Widnall  (Ref.  37)  and  of 
Platzer  et  al .  (Ref.  38)  up  to  n=5  for  both  reduced  frequencies 
(k=0  and  1.0). 

7.6.3  Salient  Features  of  BTM 

It  has  been  shown  that  the  present  method  has  the  following 
features : 
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a.  With  the  H-G  model,  the  BTM  is  computationally  efficient  and 
robust  for  unsteady  flow  computations  in  the  full  frequency 
range . 

b.  Although  the  total  number  of  panels  needed  for  BTM  is 
comparable  to  that  of  a  surface  panel  method,  the  CPU  time 
is  only  one  tenth  of  the  latter. 

c.  The  BTM  places  the  control  points  on  the  exact  body  surface 
whereas  a  surface  panel  method  usually  places  the  control 
points  on  the  approximate  surface. 
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CHAPTER  8 


Body-Wing  Combinations 
8.1  Unsteady  Lifting  Surface  Method 

Previously,  an  effective  panel  method  has  been  fully 
developed  for  oscillating  lifting  surfaces  in  supersonic  flow  HGM 
(Ref.  4).  This  analysis  has  been  further  extended  from  the 
velocity  potential  treatment  to  that  of  an  acceleration  potential 
one  so  that  the  regions  of  wake  can  be  completely  avoided  in  the 
calculation.  For  this  reason,  the  acceleration  potential  version 
of  the  HGM  (called  AHGM)  is  most  advantageous  in  coping  with  the 
unsteady  problems  for  wing-body  combinations,  in  that  only  the 
panels  lying  on  the  wing-body  surfaces  are  needed  to  be  accounted 
f  or . 

According  to  AHGM,  the  normal  velocity  and  the  unknown 
linearized  pressure  can  be  related  by  the  familiar  expression, 
i  .  e . 

|^-L-(x,y,z)  =  |  J  Ac p  K  ds  (8.1) 

A 


where 


K  =  e 


i /i(x-t)  d  f  cos(xR)  -ifi(x-v) 


-  f 

an  J 


R2  -  (x  -  v)2  -  p2[{  y  -  7)2  +  (z  -  «): 


Ac  =  c  l  —  c  u 

p  '-p  '-p 
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and  S  is  the  effective  wing  area  enclosed  by  the  intersections 
due  to  the  inversed  Mach  code  originated  from  the  control  point 
(  x  ,  y  ,  z  )  . 


8.2  Body-Wing  Formulation 
The  interference  between  the  body  and  the  wing  is 
essentially  provided  by  the  flow  tangency  condition.  In  the  case 
of  the  lifting  surfaces,  the  pressure-normal  velocity  relation  is 
the  result  of  the  tangency  condition,  i.e. 

J  AopKds  -  vN  (8.2) 

w 

in  which  vN  :  *i  the  normal  velocity,  or  the  downwash  on  the 
lifting  surface.  In  terms  of  the  body  influence  and  the  given 
wing  mode,  vH  represents  a  difference  between  these  two 

vn  =  Bw  -  (v„)B  (8.3) 

where  (vM)B  is  the  normal  velocity  induced  by  the  presence  of  the 

body;  as  related  to  Eq.  (7.1),  (vM)B  •  The  term  Bw 

represents  the  wing  mode.  Thus,  Eq.  (8.2)  becomes 

(J  Acp  Kds )  w  +  =  Bw  ,  on  the  wing  (8.4) 

w 

Note  that  in  the  RHS  of  Eq.  (8.4),  kernel  K  and  (—-*-)  are 

on  w 

evaluated  on  the  lifting  surfaces.  Similarly,  the  relationship 
on  the  body  can  be  derived,  i.e. 


,  on  the  body 


(8.5) 


(|^)3  +  (f  Ac  p  Kds  )  s  =  B 
w 

This  time,  the  kernel  K  and  ( *■)  B  are  evaluated  on  the  body 
surfaces  . 

With  Eqs.  (8.4)  and  (8.5),  the  body-wing  tangency  condition 
can  be  expressed  in  a  matrix  form, 


r  (3£) 

1 dn  '  8  8 

(— )  1 

1  an  '  8w 

{*■  } 

- 

(8.6) 

L  (K)w8 

(K)ww  J 

'Acp» 

B. 

dE 

where  the  kernel  function  (-= — )  is  related  to  the  least  square 

on 

matrix  by  the  following  equation 


0E 


L  0^]  =  {  [  u  1  ( nx  }  +  [  v  ]  { n  r  }  +  [w]{ni?}}.[LS] 


(8.7) 


nx  ,  nr  and  n^  are  the  directional  consines  of  the  body  panels; 

Bb  and  Bw  are  the  given  body  and  wing  modes  respectively.  Eq. 
(8.6)  is  solved  to  yield  the  Acp  value  on  the  wing.  The 
velocities  on  the  body  can  be  evaluated  at  the  control  points 
from  Eqs.  (7.13,  7.14  and  7.15).  Consequently,  the  pressure  on 
the  body  can  be  determined  from  Eq.  (6.18). 

8.3  Pressure  Coefficients 

8.3.1 _ Unsteady  Pressure  Coefficient  for  the  Body 

Eq.  (6.18)  provides  the  explicit  unsteady  pressure  formula 
in  tne  cylindrical  coordinate  for  an  oscillating  elastic  body. 
For  body-wing  configurations,  it  is  convenient  to  rewrite  this 
expression  in  the  Cartisien  coordinate  on  the  body  surface. 
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Cp  =  -2  S0-{(l  +  #0x)(^lx-2zBg"#0x)  (8.8) 

+  ^oy^ly  +  ^02^12 

+  g’*o2 

+  ik(  4>l  +  0O  z  g  -  zB  g’  x  )  } 

where  zB  is  the  control  point  location  of  the  body  panel  and  g  is 
the  given  body  mode. 

8.3.2 _ Unsteady  Pressure  Coefficient  for  the  Wing 

Eq.  (8.8)  can  be  applied  to  account  for  the  pressure  on  the 
lifting  surface  due  to  the  body-wing  interference,  i.e. 

ACp  =  -2So«{*Xx(l+*Ox)+#ly#Oy-*lz*0l  +  iM1)  (8.9) 

The  steady  velocities  <60X,  y  ,  ^o:  anc*  tbe  So  term  are  to  be 

evaluated  according  to  the  effective  body  panels  within  the 
inversed  Mach  cone.  Eq.  (8.9)  can  be  further  simplified  as 

Acp  =  -2S0-Olx  +  ik#  L  }  (8.10) 

since  the  second  order  terms  in  the  bracket  can  be  ignored. 
Further  simplification  of  Eq.  (8.10)  results  in  the  linearized 
pressure  formula 

icp  =  -2  {0lx  +  ik0L}  (8.11) 

which  is  essentially  the  pressure  formula  of  the  thin  wing. 

In  passing,  we  note  that  to  include  the  factor  S0  might 
upgrade  the  accuracy  in  the  higher  Mach  number  range  as  it  was 
previously  noted  by  Van  Dyke  (Ref.  10). 
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8.4  Generalized  Forces  and  Stability  Derivatives 


The  generalized  forces  Q;J  can  be  expressed  in  terras  of  the 

1 

uns :eady  pressure  coefficient  cp  and  of  the  mode  shape  g  as, 

1  f  f  i  <  J  ) 

Q  i  j  -  g -  J  J  CP  (nxg/(I>z  -  g(I)n2)ds 

S 


where  S  is  the  wing  plus  the  body  surface  area;  nx  and  n2  are  the 
normal  components  in  the  x  and  z  direction  to  the  wing  and  body 
surfaces,  respectively;  z  is  the  coordinate  of  the  control  point; 

l  <  j  ) 

Cp  is  the  unsteady  pressure  coefficient  due  to  the  J-th  mode; 

g( 1 >  is  the  I-th  mode  of  oscillation;  and  Sref  represents  the 

maximum  area  or  the  based  area  for  the  body-alone  cases,  whereas 

it  represents  the  wing  planform  area  for  the  wing-body 

combination  cases. 

The  stability  derivatives  are  related  to  the  generalized 
forces  by  the  following  formulae: 


CLa  =  Re (Q! 2 ) 


+  c 


N  q 


LAI 

k 


+  c 


-  LJQ 


m  q 


22 


kL 


where  L  is  the  body  length  for  body  alone  cases,  and  the  wing 
mean  chord  for  body-wing  combination  cases;  k  in  this  instance  is 
the  reduced  frequency  based  on  L.  Qla  and  Q22  are  the  12  and  22 
components  of  the  generalize  J  f  nee  matrix  representing  two 


rigid-body  modes.  The  former  is  the  plunging  mode  and  the  latter 
the  pitching  mode. 


8.5  Results  and  Discussion 

To  verify  the  present  results  for  wing-body  combinations  in 
the  low-frequency  limit,  measured  stability  derivatives  in  Refs. 
44  and  45  are  selected  for  comparison.  Two  configurations  are 
selected  as  shown  in  Fig.  58: 

Case  A:  Aspect  Ratio  AR  =  3.C,  Triangular  Wing-Body 

Case  B:  Aspect  Ratio  AR  -  3.0,  Swept  Wing-Body 

To  compute  the  stability  derivatives  for  these 
configurations,  a  bundle  of  20  triplet  lines  is  used  for  the 
body.  With  10  segments  along  the  body  axis,  this  amounts  to  200 
control  points  on  the  body  surface.  However,  only  50  panels  are 
required  for  modeling  the  wings. 

Figs.  59  and  60  show  the  computed  static  and  damping  moments 
for  the  wing-body  configurations  in  Cases  A  and  B  respectively. 
Good  agreements  are  found  between  the  present  results  and  the 
measured  data  in  the  overall  Mach  number  ranges  at  two  pitching- 
axis  locations,  xG  =  .25c  and  .35c  (c  is  the  chord  length).  For 
damping  moment  calculations,  one  notices  that  the  predicted 
results  of  the  analytical  theory  (Ref.  46)  deviate  further  from 
the  measured  data  than  the  present  results.  This  is  due  to  the 
fact  that  an  analytical  solution  of  a  rectangular  wing  is  used  to 
approximate  that  of  a  delta  and  a  swept  wing  planforms. 
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CHAPTER  9 


CONCLUSION 

Several  unsteady  supersonic  methods  have  been  developed  for 
computation  of  unsteady  aerodynamics  for  elastic  bodies  of 
revolution,  asymmetric  bodies  and  body-wing  configurations. 

These  methods  include  the  Harmonic  Potential  Panel  (HPP)  method, 
the  Bundle  Triplet  Method  (BTM)  and  the  combined  method  of  BTM 
and  Harmonic  Gradient  Method  (HGM)  for  body-wing  combinations. 

All  methods  are  based  on  the  Harmonic-Gradient  ( H-G )  model  in 
order  to  obtain  accurate  solutions  in  the  full  frequency  domain 
and  the  lower  Mach  number  range.  Specific  conclusions  of  the 
present  development  can  be  summarized  as  follows: 

1.  For  bodies  of  revolution  in  oscillation  the  proper 
choice  of  the  coordinate  system  has  been  subject  to  some 
controversy  in  the  past.  The  formulation  in  the  wind-fixed, 
body-fixed,  and  the  pseudo-wind-fixed  systems  are  presented.  The 
solutions  of  the  linearized  equation  in  the  wind-fixed  system, 
and  of  the  first-order  equation  in  the  body-fixed  and  the  pseudo- 
wind-fixed  systems  are  compared  for  a  cone  in  rigid  and  bending 
oscillations.  It  is  shown  that  when  the  wind-fixed  system  is 
used  in  a  straight  forward  manner  a  singularity  appears  at  the 
body  apex  and  contaminates  the  solution  downstream. 

2.  Further  investigation  of  this  problem  shows  that  such 
apex  singularities  also  appear  in  the  body-fixed  system  if  the 
linearized  equation  instead  of  the  first  order  is  solved.  Also 
it  is  shown  that  the  origin  is  due  to  the  use  of  cylindrical 
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coordinates.  When  the  problem  is  correctly  formulated  in  a 
conical  coordinate  system  the  pressure  coefficient  is  shown  to 
remain  regular  at  the  body  apex  (see  Appendix  C). 

3.  The  computational  procedure  for  the  mean  flow  has  been 
developed  for  both  the  linear  and  nonlinear  potential  equations. 
Unlike  the  unsteady  computation  procedure  for  wing  planforms,  the 
steady  mean  flow  solution  enters  into  the  boundary  conditions  and 
pressure  coefficient  for  the  unsteady  flow.  Hence,  it  is 
important  to  investigate  the  influence  of  the  steady  mean  flow 
upon  the  unsteady  pressures  as  well  as  the  flutter  boundaries. 

We  believe  that  it  is  the  first  time  that  nonlinear  effects  have 
been  included  in  an  unsteady  panel  method.  Computed  steady 
pressures  and  stability  derivatives  using  the  present  nonlinear 
method  are  found  to  agree  well  with  the  known  results.  Unsteady 
nonlinear  results  show  that  the  influence  of  the  mean  flow 
nonlinearity  is  important.  For  a  slender  body,  nonlinear  effects 
become  more  apparent  in  the  high  Mach  number  range.  Our  unsteady 
nonlinear  results  are  in  better  correlation  with  measured  data 
and  unsteady  exact  theories  than  all  other  results  using 
different  methods. 

4.  For  a  cone  in  plunging  and  pitching  oscillations  the 
flutter  boundaries,  as  determined  by  the  linear  and  nonlinear 
methods,  are  compared  to  those  obtained  by  experiment  and  other 
theories  in  Figs.  33  and  34.  With  reference  to  the  measured 
data,  all  of  the  computed  flutter  boundaries  are  conservative 
predictions.  The  computed  boundaries  become  more  conservative  in 
the  order  of  present  nonlinear  method,  present  linear  method  and 


slender  body  theory.  By  contrast,  all  results  as  predicted  by 
previous  theories  failed  to  provide  a  conservative  flutter 
boundary . 

5.  The  aerodynamic  damping  coefficient  of  a  cone-cylinder 
in  the  first  two  free-free  bending  modes  of  oscillation  and  of 
the  Saturn  SA-1  model  in  the  first  free-free  bending  mode  of 
oscillation  are  compared  to  those  obtained  by  experiment  and 
other  theories  in  Figs.  29  and  31.  For  both  configurations  all 
the  computed  results  establish  close  trends  in  the  measured  data. 
By  contrast,  all  results  as  predicted  by  previous  theories  are  in 
considerable  discrepancy.  For  both  configurations  the 
aerodynamic  damping  is  found  positive,  and  therefore  it  provides 
stabilizing  effects  within  the  Mach  number  range  considered. 

6.  With  the  H-G  model,  the  BTM  is  computationally 
efficient  and  robust  for  unsteady  flow  computations  in  the  full 
frequency  range. 

7.  Although  the  total  number  of  panels  needed  for  BTM  is 
comparable  to  that  of  a  surface  panel  method,  the  CPU  time  is 
only  one  tenth  of  the  latter. 

8.  The  BTM  places  the  control  points  on  the  exact  body 
surface  whereas  a  surface  panel  method  usually  places  the  control 
points  on  the  approximate  surface. 

9.  To  validate  the  present  method,  numerical  examples  have 
been  studied  for  various  three-dimensional  configurations.  These 
include:  steady  pressures  on  three  asymmetric  conical  bodies, 
generalized  forces  for  cylinderical  panel  flutter  and  stability 
derivatives  for  bodies  and  wing-body  combinations.  It  is  found 
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that  all  computed  results  are  in  good  agreement  with  existing 
theoretical  results  and  measured  data. 

10.  The  combined  use  of  the  developed  AHGM  and  the  BTM 
allows  the  wing-body  surface  to  be  the  only  computed  domain. 
Moreover,  the  body-wing  formulation  can  directly  yield  the 
unsteady  pressure  on  wing  planforras.  Hence,  the  present 
computation  procedure  clearly  leads  to  the  development  of  a  user- 
friendly  program  code. 

Finally,  we  conclude  that  all  methods  proposed  have  been 
thoroughly  validated  with  existing  theories  or  measured  data.  It 
is  believed  that  all  methods  developed  thus  far  can  become 
effective  tools  for  supersonic  aeroelastic  applications  to 
realistic  configurations. 

Throughout  three  phases  of  this  development,  most  of  our 
work  has  been  presented  in  the  A I AA  meetings  as  well  as  appeared 
in  various  kinds  of  AIAA  Journals,  (See  Refs.  47  to  54).  For 
technology  transfer,  one  report  along  with  the  developed  program 
code  has  been  forwarded  to  Army  Missile  Command  (MICOM)  in 
Huntsville,  Alabama  in  1986  (Ref.  55). 


APPENDIX  A 


LINEARIZED  EQUATIONS  IN  THE  BODY-FIXED  COORDINATE  SYSTEM 
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appendix  a 


The  body-fixed  coordinate  system  can  be  related  to  the  wind- 
fixed  coordinate  system  of  Fig.  3  by  the  following  transformation 

x  =  x’  -  |^r(x’  ,  t’  )  z’  (A.  1) 

y  =  y’ 

z  =  z’  +  h(x* , t’  ) 


where  h(x' , t  ’  )  is  the  instantaneous  normal  displacement  of  the 
body  and  only  linear  terms  in  h  are  retained  since  it  is  assumed 
to  be  smal 1 . 

The  operators  7,  72  and  7*  of  Eqs.  (2.2)  and  (2.3)  in 
curvilinear  coordinates  are  given  by  (see  Ref.  39) 


1_  d»_ 

hi  a9i 


ai 


1 _  3y 

h2  aq2 


,  l  a» 

h3  aq3 


(A.  2) 


7* A  = 


^1  ^2  ^3 


aqi(h2h3Al)  +  3q2  ^  h3  A2  ^ 


dq  (hl  h2  "3 


A,  ) 


(A. 3) 


and 


72  if  = 


hi h2h3 La<li 


.(  jla 


a* 


ht  aqz 


dq. 


,h,h,  dt 

h2  aq: 


d  q 


-(  hlh2. 


h.. 


a*  . 

dq3 

(A. 4) 


where  qx  ,  q2  and  q3  are  the  curvilinear  coordinates,  ax  ,  a2 ,  and 
a3  are  unit  vectors  along  these  coordinates,  hx  ,  h2  and  h3  are 

the  matrix  coefficients  defined  as 


ax’2  av’  2  az’2 

h2  =  (!£-)  +  (^-)  +  (|§-) 


dq, 


dq, 


aq, 


(A. 5) 
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and  x’,  y ’  and  z  ’  are  the  cartesian  coordinates. 


In  the  present  case 


hi  =  1 
h,  =  1 


h3  =  1 


<3  2  h 

dx2 


(  A.  6) 


where  the  nonlinear  terms  in  h(x,t)  have  been  neglected. 

The  total  potential  0  in  the  body-fixed  coordinate  system 
can  then  be  expressed  as 

Q(x,y,z,t)  =  x  +  hx  (  x  ,  t )  2  +  <t(x,y,z,t)  (A. 7) 

where  4>  (  x  ,  y  ,  z  ,  t )  is  the  perturbation  potential.  From  Eqs.  (A. 2), 
(A. 4)  and  (A. 6)  and  retaining  only  linear  terms  in  h  we  obtain 

7Q  =  (1  +  *x  -  hxxz*x)  ex  +  iyey  +  (hx  +  *z)e2  (A. 8) 

Q  -  hvvvz  +  v  -  2h_  z$„  „  -  h„„„z  -  h„„„z$„ 

"AA  XX  XX  XX  XXX  XXX  X 

+  +  hxx$z  +  Qzz  (A. 9) 


The  time  derivative  in  Eqs.  (2.2)  and  (2.3)  in  the  body- 
fixed  system  is  given  by 


<3  d  dx  d  dz  d 

at  ’  It  +  aF  dx  +  dT7'  dz 

d  ,  a  ,  a 
at  x  t  2  0X  +  t  dz 


(A.  10) 
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Therefore 


ao 


at  ’ 

a2  q 


=  h, 


<t.  -  h  x  ,  z  (  1  +  $  )  +  h .  $  _ 


(  A  .  1 1  ) 


dt‘—  =  h  x  .  ,  z  *  * .  .  -  h;< ,  t  z  (  1  +  ♦  x  )  -  h  :<  t  z  <t  x  * 

+  h , .  ,  <J>  z  +  h ,  $  2  t  A  .  12} 

We  can  now  compute  all  the  terms  of  Eqs.  '2.2)  and  ,2.3 
the  body-fixed  coordinate  system.  If  only  linear  terms  in  <t> 
retained  the  following  linear  equation  is  obtained 

*  „ 


2  h  z  $ 

4-11 XX  ^  *X  X 


h  z 

X  XXX 


h 

y  y 

X  X 

„ 

2  2  2 

z*xt  +  ht 

t*. 

+  h,*Jt 

X  x  Z*X  t  + 

2hx 

t*,  +  2hx« 

2  t 

-  2hxtz$xx  +  2ht#X2  +  ♦xx(l~2hxxz)  h  x  x  x  z x 

+  hXx*2  +  2hx*xZ  ]  (A.  13) 


Now,  letting  h(x,t)  =  g(x)<5(t)  =  J0e'  kt  g(x),  and 

♦  (  x  ,  y  ,  z  ,  t )  =  <p0  (  x  ,  y  ,  z  )  +  5o  e  1  *  f-  4>l  (  x  ,  y  ,  z  )  ,  with  50«1,  and 

substituting  these  into  Eq.  (A. 13)  yields  the  following 
equat ions: 


{  1  -M2  ;  <t>. 


o  y  y 


0 


i A . 14) 


for  order  of  one  and, 


)  l  n 
are 
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*lxx  -  ^lyy  +  tlzz  ~  2ikMl*lx  + 


-  ^[k^'zfi:<  -  k  2  g  0 ,,  ,  -  2ikg’’z0Ox  -  2g’ik?> 


0  2 


g 


2g’  ikzflu*  2gik0Oxz  -  g’’’z^0x  -  2g’’z?>0xx 
2g’^oxz]  "  g”*o2  +  2g”^oxx  +  g”’^ox  f  A.  15) 


for  order  of  5 .. 


In  cylindrical  coordinates  Eqs.  (A. 14)  and  (A. 15)  become 


(1-M2)  0 


0  x  x  r  r 


+  — -  <p 

i  r  r  2 


0  ee 


A.  16) 


'  1-M2  )  +  d> ,  +  -  <i ,  -  - —  6  ,  -  2  ikM2  <b  ,  +  k2M2tf 

'  ^Ixx  '1  rr  r  "lr  r2  ^  1  qq  o«v  1  x  A  “oo^l 

=  M^[(  k2  g’  r?>0  x  -  k2  g^0  r  -  2ikg’’r0o>  +  2g’ik0Or 

-  2g’ikr*0xx  +  2gik^0  x r  -  g’”r*0x  -  2g’’r*0xx 
+  g’’0or  +  2«’^oxr)  cos0  +  ( k2^O0  -  2g’ik^O0 

-  2  g  i  k  x  Q  -  g”^o9  -  2g’^OX0)  p  sine] 

-  g’  ’  lfJrcose  -  -  0O0  sin©)  +  2g  ’  ’  r?>0  x  x  cose 

+  2  g ’  ’  ’ r^0xcose  (A.  17) 


For  a  body  of  revolution  the  mean  flow  potential  <f>0  is 
independent  of  6,  and  the  cross  flow  potential  can  be  expressed 
as  *  j  ( x , r , e)  =  ^  fx,r)cose.  Thus  Eqs.  (A.  16)  and  (A.  17)  become 

1 


■1-M*)  *oxx  - 


0  r  r 


-  <t^  =  0 

^  0  r 


(  a  .  1 8 ; 
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(l~^V  <>lxx  +  *lrr  +  7  ^lr  ~  72"  *1  ~  2ikMl*ix  + 

-  M2[k2g*r^0x  -  k2g*0r  -  2ikg”r*0x  +  2g’ik* 

-  2g’ikr^0xx  +  2gik^0xr  -  g’”r*0x  -  2g”r*0 

+  *"*or  +  2S^Oxrl 

+  2g”^0xx  +  g”’^Ox  -  g’^or 


k2M^x 

0  r 

X  X 


(  A  .  1  9  ) 
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APPENDIX  B 


EVALUATION  OF  THE  UNSTEADY  POTENTIAL  AND  VELOCITIES  APPLYING  THE 
•  HARMONIC  GRADIENT  MODEL 
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APPENDIX  B 


The  integral  solution  of  the  homogeneous  unsteady  wave 
equation,  i.e.  Eq.  (2.8)  reads 

1  rx~Pr  a 

M*»r)  =  "27  J  F(f)  fp  K(x-S,/Jr)d?  (B.l) 


where  the  kernel  function  K,  as  defined  in  Eq.  (4.2)  is  given  by 


K(x-<  ,  /Sr)  -  e 


i/i(x-S  )  cos  XR 
_ 


(B.  2) 


where  R  =  /( x- S ) 2  -pz  r2  ,  n  =  kM^/£2  ,  X  =  kM^//®2 

and  F(5)  is  the  dipole  strength.  The  value  of  F(f)  at  f=0  is 
zero  in  order  to  satisfy  the  Mach  wave  condition. 

For  integration  of  Eq.  (B.l)  it  is  convenient  to  use  the 
relative  distance  x-S  as  the  independent  variable  instead  of  t. 
Hence,  we  define  x0  =  x  -  S  and  df  =  ~  dxQ 
where  at  5  =  0  ,  x0  =  x 

and  at  (  =  x  -  fir  ,  x0  -  p r 

and  Eq.  (B.l)  becomes 


(x, r) 


i  rfr 

=  27  F<x-*o> 

*  V 


,"i^xo  a  , cos xrq 

<3r '  Rn  '  ' 


(B.3) 


where  R0  =  /x20 -/J2  r2  . 

Integrating  Eq.  (B.3)  by  parts  and  applying  the  Mach  wave 
condition  at  x-p r=0,  we  obtain 


*i<x.r)  =  F/x_xo  )  e  1|UX°  IpS  (x0  .  Ar) 


where  S ( x 


f 


o  .0**)  =  J 


fir 


J^-[F(x-x0)e  1MX°]  |^S(x0  ,  (5r)dx0J 


(  B  .  4  ) 


cos  x/t2  -ft2  r2 


-d  t 


ftr  / r2  -ft2  r2 


The  first  term  of  Eq.  (B.4)  is  zero  at  the  upper  and  the 
lower  limits  since  F(0)=0  for  a  pointed  body  and  S(ftr,ftr)=0  along 
the  apex  Mach  wave.  Consequently,  0L(x,r)  can  be  expressed  as 


<t>i  ( x 


-  -h  J 


fir 


dx, 


-[  F(x-x0  ) 


i^xc 


It® 


(x0,ftr)dx0  ( B . 5  ) 


As  discussed  in  Section  4.2,  to  solve  for  Eq.  (B.5)  the 
dipole  strength  function  is  discretized  into  elements  along  the 
body  axis  and  the  potential  is  evaluated  at  the  discrete  field 
points  (  Xj  ,  Tj  )  . 


*i (*j «rj > 


[ F,  (Xj~x0)e 


"iMX0  j 


S(x0 ,ftrv ) dx0 


(  B  .  6  ) 


where  !1  =0  and  !j  +  x  =Xj  -ftr^  . 

Introducing  the  Harmonic  Gradient  model  (as  discussed  in 
lection  4.2), 


1 _ r  p  (x  -x  )  e  i/zX< 

ax0L(‘uJ  x°; 
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Eq.  (B.7)  amounts  to  stating  that  the  dipole  strength  is 


+  — £-)  ( e  1//(XJ  X°}-  1)  -  ?->-(x  -x0  )  (B.8) 

Ifl  XfL  J  u 


In  the  limit  of  k-*0  (  or  Eq.  (B.8)  is  reduced  to 


F, 


-(Xj-Xq)2  -  b!  (Xj-Xq) 


(  B  .  9  ) 


Now,  the  doublet  strength  model  as  given  by  Eq.  (B.9)  is  the 
same  as  the  one  used  by  Tsien39  for  the  steady  angle  of  attack 


case . 


Substituting  Eq.  (B.7)  into  (B.6)  yields 


1  f  rxj  *i  +  i 

(xj  >  rj  )  =  "27  l  J  _  t  a>  ( xj  -xo 


)  +  bj  ]  e 


-inx. 


i=i  xJ"ei 


af-  S(x0  ,^frJ  )dx0 


(B.9) 


In  order  to  determine  S(x0,/9rj)  it  is  convenient  to  apply 
the  following  transformation: 

Let  t -/Sr  ■  cosha  ,  and  hence  dr  =  /5r ^  sinhtrdcr 


where  at  , 


and  r=x 


o  > 


cr-0 


<r  =  cosh-  1 


a  a  fCosh" 1  x0  / 

yp-S  ixQ,pri)  =  jp-  j  cos(\/irjsinhff)d(7. 


By  applying  Leibnitz’s  rule, 
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xn  cos  (  \/x2  -pz  r,2) 

r j  rj2 

(B.  10) 

(•cosh'  1  x0  //Sr. 

\/3  sinhtr  sin(X/8r  sinh<7)d<7 

J  o 


Using  Lashka’s27  exponential  series  substitution,  namely 


/1  +  u2 


H 

l  ane'nCU 

n  =  1 


where  the  constant  "c"  and  the  coefficients  "an"  can  be  found  in 
Ref.  27.  Letting  x0  =£r_j  cosht ,  and  integrating  Eq.  (B.10)  term  by 
term  yields 
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dr 


-S  (  r  ,  p  r  J 


cosht  cos  ( \pr  ,  s  inhr  )  cos(Xjgr,sinhT) 

x^r  j  P 


r  s  inhr 


N  -ncsinhi ,  „ 

_ _ _ \fi 


l 


(  \p r  ,  ) 2  +n2  c2 


n  =  l 


[  -  n  c  sinCX/JTjSinht)  -  X/SrJcos(X(5rjsinhr)] 


x£r 


_  V  an  (-^r,  )  M 
L  (X^rj)2+n2c2 


(B.  11) 


n  -  1 


By  substitution  of  Eq.  (B.ll)  into  Eq.  (B.9)  we  obtain 


*i (Xj - rj )  =  "27 


cosh-1  (  — -1-— 

/*rj 


‘  =  1  J  cosh-  1  -*■■<-) 

j 


[at  (Xj-^TjCoshr)  +  bt  ] 


^-i/z^r  j  coshi  f-/ff  coshr  cos  (  X/STj  s  inhr  ) 


+  sinhr  p 


-1  +  cos  (  X/STj  s  inhr  ) 


i  -ncsinhr 

I  (Mr,  ) 2  -.  {  nc )  5  I-"  c  »in(MrjSi»hr) 


n  =  1 


r  a  (  -  x2  6Z  r2  )  ^ 

-x^rj  cos(X/ffrJ  sinhr)  ]/»rj  X  -  £  (x/8rJ  )2+(nc)2J/d 


n  =  1 


(B. 12) 


Since  all  the  terms  of  the  integrand  in  Eq.  (B.12)  are 
regular,  they  can  be  integrated  numerically  by  using  a  Gaussian 
quadrature  formula. 

The  velocities  in  the  x  and  r  directions  are  determined  by 
taking  the  derivative  with  respect  to  x  and  r  in  Eq.  (B.3). 
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After  the  derivatives  are  taken  the  dipole  strength  function  is 


discretized  and  the  Harmonic  Gradient  model  applied. 

The  x-component  of  the  velocity  can  be  expressed 


a  .  ,  .  1 

air*.  'ri>  =  - 


cosh  1  ( * 1  — 1 ^ ) 
1  J  cosh'1  — ^-L-) 


-  iu.Br,  coshr 
a,  e  J 


^-/Jcoshr  cos  (  V^r  j  s  i nhr  )  +  /Ssinhr 


1 


[  -nc 
u 

~  I 

n  = 


+  cos  (x/Sr^s  inhr  ) 


sin ( sinhr ) 

an  (~\2  p2  r ,  2  ) 

(  *P )  2  +  (  nc  )  2_ 


} 


-ncsinht 

an  g _ 

(  X/5r  j  j2  +  tnc)2 

n  =  1 

cos(  X/STj  sinhr)  ]  pr  ^  X 

dr  ( B . 13) 


H 

l 


Again  Eq.  (B.13)  can  be  integrated  numerically  by  using 
Gaussian  quadrature  since  all  the  terms  are  regular. 

The  r  component  of  the  velocity  becomes 


d  i  f  rxj  (i  +  i 

a7~^i  (*j  .r,  )  =  -  27  1  J  _  fai(xj-xo 

i  =  1  X  j  ~  *  i 


)  +  b,  ]e 


-1AXC 


afT  S(x0  ,pr3  )dx0 


(B. 14) 


next,  letting  x0  =  /JTj  cosht  and  applying  the 
differentiation  with  respect  to  Tj ,  to  Eq.  (B.10)  yields 
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I > 


d _  cosht  cos(  X|8r,  sinhi  ) 

3rj  TjSinht 


x^sinhcr  s i  n  ( X/5r  s  i nha )  da 


cosht  cos  ( \/8r  ,  sinht ) 
r^sinht 

cosht  \/3  s  i  n  (  \/3r  ,  s  i  nhr  ) 


X2 /?2  s  inh2  crs  in  (  x£r  s  inho-)  do- 
Jo 


(B.  15) 


When  the  transformation  sinho=u  is  introduced,  the  integral 
term  of  Eq.  (B.15)  becomes 


•sinht  2 

l  X2  pz  —■  s  in  (  X/Jr  .  u  )  du 

lo  /1  +  u2 


Again,  when  Lashka’s  exponential  series  substitution  is  used  this 
integral  becomes 


,  ,  /Sr, sinht  •  u  \  .  cos  (  X/8r  ,  s  inht )  ,  _ 

x2PZ  -Tt>rJ  )2  sin(X^rjSxnht)  +  - a - -  x2/*2 

5.  -ncsinht  r 

-  x2/J2  I  ) 2  +  ( nc ) 2  |_MinhT(~n  c  cos(x/JrjSinht)) 

X2/?2 


n  =  1 


+  X/Sr  j  s  in  (  X/Sr  s  inht ) 


(  nc  )  2  +  (  X^Tj  ): 


( (nc)2-( X£r  ^  )2  ) 


cos  (  X/Srj  s  i  nht )  -  2  n  c  x^rj  s  in  (  x^r^  s  inht ) 
x2/?2 


x2  £2  rj 


V  x2/S2a„((nc)2-(xffr,)2) 
Z  (  (nc)2  +  (  \firi  )2  )2 


(B. 16) 


n  =  1 


a2s 


With  Eq.  (B.16),  ^ - j-  of  Eq.  (B.15)  can  now  be  expressed  in 

°  r  j 
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a  series  form,  i.e., 


32 


S(r  nr  ^  -  coshr  cos  (  Xftr  .  sinhr  )  ,  \fi  coshr  s  i  n  (  \/S  r  ,  s  i  nhr  ) 
<3r,2  r , 2  s  inhr  7 - - 


^-SI?hT  sin(X|5r  sinh:)  -  ^.,s  (  Mr  ,  s  inhr  ) 

i  J  r*  2 


\2  p2  y  ?n  e 

L  {\fir 


-ncs i nhr 


) 2  +  ( nc ) 2 


n  =  1 


j^c  inh  r  (  -n  c  cos  (  \/3  r :  s  i  nhr  )  +  x/?r s  i  n  (  s  i  n hr  ) 

\2  RZ 

+  I  nc)2  \/3Vj  )2  [  (  (  nc)2  -  (  \  ^Tj  )2  )cos(  X^rj  sinhr  ) 

-  2  n  c  X/Srj  s  in  (  XjSrj  s  i  nhr  )  ]  -  ply 

-  X  x2^an((nc)2-(x/?r,)2) 

Z  (  (nc)2  +  (  X/Jr,  )2  )2  (B.17) 


n  =  l 


Finally  where  Eq.  (B.17)  is  substituted  into  Eq.  (B.14)  and 
after  some  arrangement,  we  obtain 


cosh-1  (  —■ >■ . *  >  t  x  ) 


1  =  1  J  cosh'1  (1 — !->-) 
Pri 


[  a  j  (Xj-jJrjCoshtJ+b, 


-i/ij8r  coshr  .fcoshr  . 

e  J  /S | - — — costx^r^  sinht) 

+  \p  coshr  sinhtsin(  X(8r0  sinhr) 


\p  sinh2r  s  i  n  (  X/S  r  s  i  nhr  )  -  —  T. 

J  r  , 


cos  (  y  fi  r  j  sinhr  )  +  V2j«2rjsinhT 


H 

l 


a^e 


-ncsinhr 


(Mr,  )  2  +  (  nc) 


n  =  X 


■j-  (sinhr(-n  c  cos  (  X/S  r  j  s  inhr  ) 
x2  ($2  r  ,  s  inht 


.  ^sxnf^sinh,)]  -  {^XqT-) 

[  ((nc)2-(  Mrj  )2)cos(X/Jrjsinhr) 

-  2  n  c  s  in  (  X/Jrj  s  inhr  )  ]  |  dr 


(B. 18) 


All  the  terms  in  E  q . 
integrated  numerically  by 


(B.18)  are  regular  and  can  be 
a  Gaussian  quadrature. 
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APPENDIX  C 


ON  THE  APEX  SINGULARITY  FOR  OSCILLATING  POINTED  BODIES:  WIND 
FIXED  VERSUS  BODY-FIXED  COORDINATE  SYSTEMS 
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APPENDIX  C 


It  has  been  known  for  some  time  (see  Refs.  17,  41,  42)  that, 
in  the  cross  flow  potential  for  an  oscillating  body  in  the  wind- 
fixed  coordinates  system,  a  singularity  at  the  apex  arises.  The 
origin  of  this  singularity  is  the  second-order  derivatives  of  the 
mean  flow  that  appear  in  the  boundary  conditions  and  unsteady 
pressure  coefficient  when  these  are  evaluated  at  the  mean 
position  of  the  body  surface.  For  a  body  in  rigid  oscillation 
with  respect  to  x ,,  the  tangency  condition  and  the  unsteady 
pressure  coefficient  are  respectively  given  by 

*lr  ~  H’Ci,  =  1  *  +  (X-XQ  )  (*0  r  r  "  R  ’  <P0  x  r  “  ik) 

+  (x-xG)(^0rr  -  R’^0xr  ~  ik)  at  r=R(x)  (C.l) 


1  -  2S0[(1  +  *0x)*lx  +  t0rtlr  +  ik*l 

-  (  X-XG  )  (  <pQ  x  r  (  1  +  *0x)  +  *0rr*or)] 


at  r  =  R(x)  ( C . 2  ) 


where  S0  =  [1  -  ^  M^(2*0x  +  *2X  +  #§x)]1/r  1,  k  is  the 

reduced  frequency,  <pQ  the  mean  flow  potential  and  <p  1  the  unsteady 
flow  potential.  The  second-order  derivatives  of  <p0  can  be  better 
shown  to  be  singular  from  slender  body  theory.  The  velocity  <pQ  r 
in  the  slender  body  theory  (Ref.  5)  ic 
RR  * 

given  by  <p  0  r  =  — — ,  and  its  derivatives  with  respect  to  x 


and  r  by 

R ’ 2  +  RR  ’  ’  ,  ^  _  RR’ 

r  x  “  “  ana  '  o  r  —  ^ 


(  C  .  3  ) 
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On  the  body  apex 


these  terms  behave  like 


il  V  V  } 

Rlol  ; 


therefore 


for  a  nonvanishing  R  ’  (  0  )  they  are  singular. 

In  the  body-fixed  coordinate  system  the  singularity  at  the 
apex  has  passed  usually  unnoticed  because  the  first-order 
equation  instead  of  the  linearized  one  has  been  solved  in  this 
system  and  then  the  second- order  derivative  terms  do  not  appear 
in  the  formulation. 

For  a  body  in  rigid  oscillation  with  respect  to  x .  the 
complete  formulation  in  the  body-fixed  system  is  given  by: 
Governing  equation; 

U-MV  *lxx  +  *lrr  +  ~  *lr  -pr  -  2ikM^lx 

=  Ml^k2  )*or  >  +  2ik(^0r-r^0xx  *  (X-Xs  )^Qxr  > 

+  2^Oxr]  ( C . 4 ) 

Boundary  conditions; 


^i  -  -  0lr  =0  at  x-p r<0  I C . 5 ) 

$ i r  -  R’#ix  -  “1  -  ik ( x-xG fRR ’ )  at  r=R(x)  (C.6) 

Pressure  coefficient; 

C?  =  "  2So^lx(1  +  ^x)  +  +  i  k  (  ( X  -  X  3  )  *  3  r 

-  R0O  x  -1-  <t>x  )  ]  at  r  =  R  (  x  ;  (  C  .  7  ) 


If  the  solution  of  the  linearized  equation  in  the  body  fixed 
is  attempted  we  need  to  account  for  the  right-hand-side  terms  in 


Eq.  (C.5).  This  can  be  done  through  the  following  particular 
solution 


0,p(x,r)  =  r  0O  x  -  (x-xG)*0r  (C.8) 

Therefore  the  solution  of  Eq.  (C.5)  can  be  expressed  as 
^(x.r)  =  0j(x,r)  +  ^P(x,r)  where  is  the  solution  to  the 
homogeneous  equation. 

By  substitution  of  this  equation  into  Eq.  (C.7)  and  (C.8)  we 
obtain 


*ir  -  R’*1x  =  "1  -  ik(x-xG  +  HR* ) 

-  <f>Qx  -  R^oxr  +  (x'xG)‘(^Orr  “  8’^qu) 

+  R’  (^ou'^or)  at  r  =  R(x)  ( C .  9 ) 

and 

CP  =  _2S0  I*?*  (Wox  )  +  ^or^lr  +  ik*l  +  R^Ou 

"  (X"XG  H^on  (  1  +  *0x  )  +  *Or*Orr) 

+  R^Ox  *Oxx  +  ^Or  ^Oxr  >>  (C.10) 

It  can  be  noticed  that  when  the  linearized  equation  is 
solved  in  the  body  fixed  system  also,  second-order  derivatives  of 
the  mean  flow,  and  therefore  singular  terms,  appear  in  the 
tangency  condition  and  pressure  coefficient. 

In  terms  of  slender  body  theory  Hoffman  and  Platzer12  have 
shown  that  although  the  velocity  #lr,  in  the  wind-fixed 
coordinate  system  is  singular  at  the  body  apex,  the  pressure 
coefficient  is  regular.  In  the  body-fixed  system  Revell21  also 


has  shown  that,  within  the  approximation  of  the  slender  body 
theory,  the  pressure  coefficient  as  obtained  from  the  solution  to 
the  linearized  equation  for  the  unsteady  cross  flow,  also  remains 
regular. 

In  terms  of  no t-s o-s 1 ender  body  theory,  this  apex 
singularity  has  not  been  investigated  before.  In  this  Appendix 
the  formulation  in  the  wind-fixed  system  for  unsteady  flow  is 
presented  in  a  conical  coordinate  system  (see  Fig.  58)  and  shown 
that  in  this  conical  coordinate  system  the  apex  singularity  is 
totally  removed.  For  pointed  bodies  the  body  apex  can  be  treated 
as  a  cone  and,  therefore,  the  assumption  of  conical  flow  and  the 
use  of  a  conical  coordinate  is  justified  at  the  apex. 

The  conical  and  cylindrical  coordinates  are  related  by  the 
following  equation 

s2  =  x2  +  r2  x 

}  (C.ll) 

tan  v  =  r/x 

The  mean  flow  velocities  in  the  x  and  r  directions  are  given 
by  ( see  Ref.  43 ) 


0Ox  =  -  a  cosh-1  ( C°tv ) 
<fi0r  =  a  /cot2  v-pz 


(C.  12) 


where 


1 


co  t  <7 /co  t 2  <7-fSz  +  cosh  1  ( 


cot  ff 

~7 


and  cr  is  the  half-cone  angle. 
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It  can  be  seen  that  no  singularity  appears  in  this  equation. 
The  form  of  Eq.  (C.14)  suggests  that  the  solution  <fi1  can  be 
sought  in  the  form  of 


0!  (  S  ,  V  ) 


J  sn  f  n  (  V  ) 
n  =  0 


(c. 15) 


When  the  governing  equation  for  <fi  l  in  the  wind-fixed 
coordinate  system  is  transformed  to  the  (x,v)  coordinate  and  Eq. 
(C.15)  is  used  we  obtain 
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sn  '  2  [f  n  (  v  )  (  n  (  n- 1  )  (  s  in2  v-/52  c°s2  v  )  +  n(cos2v-/52sin2v) 


n  =  0 


+  n  -  — r-^-5 — )  +  f  ’  (  v  )  (  2s  in  vcos  v  (  n- 1 )  M2  +  cotv) 
sm2  v  n  - 


+  f  „  ’  (  v  )  (  p2  s  i  n2  v  +  cos2v)]  +  k2M^  ^  s n  f  n  (  v  ) 

n  -  0 


-  2ikM^[  ^  sn ”  1  [ f n ( v ) ncos v-f ^ ( v ) s in v ] ] 


n  =  0 


0 

(C. 16) 


When  Eq.  (C.15)  is  used  to  satisfy  the  boundary  conditions, 
Eq.  (C.14),  the  following  equations  and  boundary  conditions  for 
f0 ,  fj  and  fj  (where  j=2,3,...)  can  be  obtained  from  Eqs.  (C.16) 
and  (C.14) 


(-/S2sin2v  +  cos2v)  f ’ ’ (v)  -  2M^  sinvcosvf  ’  (  v  ) 


+  cotvf  ’  ( v )  -  sin2-~  f0(v)  =  0 


f0(v)  =  0 


f  ’  (v)  =  s.  (- 


sin3  a  +  /cot2  a- p 2 


at  v-fi 
(Mach  wave) 

+  ik)  at  v  =  <7 

( cone  surface ) 


(C.17) 


where  s,-  is  the  center  of  rotation  in  the  conical  coordinate 
system. 
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( -  p2 s i n2 v +  cos2 v ) f ’ ’ ( v )  +  cotvf’(v) 


-  (/52sin2v  +  cos2vcot2v)f1(v)  =  -  2  i  kM^s  i  n  v  fQ  (  v  ) 

f  j  (  v  )  =  0  at  v -fi 

f (  v  )  =  -  c  o  s  <7  f  1  +  acosh"  1  ( — — ) 

P 


(C.  18) 


.  . . -  -  ik] 

sin3  <r/cot2  <j-pz 


at  v-<r  J 


and  finally  a  recurrence  formulation  can  be  written  for  fj(v) 
with  j=2 , 3 , . . . 


f j  ( v ) [ J ( j~ 1 ) [ s in2 v-p2 cos2 v )  +  j ( cos2 v-p2 s in2 v ) 

+  j  -  — r-U — ]  +  f  ’  (  v  )  (  2sinvcosv  (  j-1  )M2  +  cotv) 

sin‘v  J  00 


+  f J  ’ ( v ) ( p2 3 i n2 v+cos2 v ) 

=  -  k2M2  f(v)  +  2  ikM2  ( f ( v ) ( j - 1 ) cos v 
**  j  -  2  **  j  -  i 


(C.  19) 


f  ’  (  v  )  s  i  n  v  ) 
j  -  i 


fj  (y)  =  o 
f’  (v)  =  o 


at  v  =  /i 

at  v  —  (7  ,  for  j>2 


It  can  be  seen  that  the  solution  to  Eqs.  (C.17),  (C.18)  and 

(C.19)  do  not  give  any  singularity  at  the  apex  and  therefore  the 
potential  as  given  by  (C.15)  will  be  regular  at  the  apex. 
However,  as  with  the  slender  body  theory,  the  velocity  <fi  x r  will 
be  singular  at  the  apex.  As  can  be  shown,  it  behaves  like  1/s. 


Ill 


We  need  to  investigate  whether  or  not  the  pressure  remains 
regular.  To  do  so,  let  us  consider  first  the  terms  + 

(lorr^or  in  ®*q.  (C.2),  which  can  be  expressed  at  v-a  as 


a  co t <7 


[l  -  a  cosh" 1  ( ~  a  cot<r/cot2  <r-p2~\ 


s  ina/cot2  <j- p2 


P 


and  by  substitution  of  the  value  of  a 


a  co t cr 


ssinaVcot2  <r- p 


2'  L 


1  - 


,  CO  t  <7  x 

1  i  *— > 


(cosh-1  p  +cot  17/cot2  <7-p2 


cotff/cot2  cr-p2  +cosh-  1  (  -  ° ) 


=  0 


and  Eq.  (C.2)  for  the  pressure  coefficient  in  the  (s,v) 
coordinate  can  be  expressed  as 


c  t  n  V 

C1  =  -2S0  { ( 1  +  0OX  )  (cos^j  s  -  — - — <f>  l  x  )  +  <f>0  r  (sinv^ 


s 


+  <t>x  } 

S  1  V 


(C. 20) 


If  Eq.  (C.15)  is  substituted  and  the  terms  are  rearranged, 
we  obtain 


=  ~2S0 


^  sn  -  1  (nf  n  (  v)  [  cosv  (  l  +  <pQ  x  )  +  sinv^0  r 


n  =  0 


-  f’  (v)  (sinx(l+0Ox  )  -  <p0rcosv))  (C.21) 

But  sinv(l+*0x)  -  ^Qrcosv)  =  -0Or  +  tanv(^0x+l)  -  0  at  v  =  cr 
by  the  mean  flow  tangency  condition.  Therefore 
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=  -2S0  2  sn  1 nf n ( v) ( cos V ( 1  +  0O x )  +  sinv  0Or)  (C.28) 

n  =  O 


Eq.  (C.28)  is  finite  everywhere.  Therefore  it  is  proved 
that  for  not-so-s lender  oscillating  bodies  the  singularity  at  the 
body  in  the  wind-fixed  coordinate  system  is  a  spurious  one  and  if 
the  problem  is  formulated  in  a  conical  coordinate  system,  the 
pressure  is  regular  at  the  apex.  The  same  conclusions  can  be 
obtained  if  the  previous  formulation  is  done  for  the  body-fixed 
coordinate  system. 

Since  the  singularity  in  the  conical  coordinate  system  is 
removed  this  solution  for  the  body  apex  can  be  used  as  the 
starting  solution.  For  the  rest  of  the  body  a  cylindrical 
coordinate  system  can  be  used  and  since  the  singularity  has  been 
removed  a  regular  pressure  coefficient  will  be  obtained 
everywhere  on  the  body. 
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APPENDIX  D 


EVALUATION  OF  THE  INDUCED  UNSTEADY  POTENTIAL  AND  VELOCITIES  BY  A 

DISTRIBUTION  OF  SOURCES 


APPENDIX  D 


As  shown  in  Chapter  7  the  induced  potential  by  a 
distribution  of  unsteady  sources  at  the  point  x4  ,  ri ,  <fim  is  given 
by 

x i ~£r i 

<t>l  (x(  ,  r.  ,  9m  )  =  |  Fn(«)B(xl  -  i,  (D.l) 

0 


When  the  transformation  S  ~  xi  -  pr  {  cosher  is  introduced  to  Eq. 
(D.l),  the  potential  <fi  L  (  x  t  ,  r  i  ,  <rm  )  can  be  expressed  as 

0 

0i  (*i  *  ri  *  Fm(xi  "  lSri  cosher)  BC^/JrJd^  (D.2) 


cosh- 1 


X  { 
^ri 


where  B(s*,/ffrj)  -  e  cos^°'  .  cos  ( \prx  s  inh^ )  . 

After  discretizing  along  the  x-axis  and  upon  applying  the 
Harmonic  Gradient  model  Eq.  (D.2)  can  be  expressed  as 


1  er 


i*  1 


(Xj  ,  r4  ,  erm  )=-£  a  j  j  [e  1MXi~e  1  r  ‘  COsha  ]  ♦  cos  (  x^r  ;  s  inher )  der 

j  =  l  (D.4) 

where  er  ,  =  cosh'  1  — — *-)  and  cr  =  cosh'1  (  X  1  - g-j-t.^-)  . 

J  P  r  i  J  *  1  pr. 


Since  the  integrand  in  Eq.  (D.4)  is  regular  the  integration 
can  be  done  numerically  by  a  Gaussian  quadrature. 

The  velocities  in  the  x  and  r  direction  can  be  obtained  by 
taking  the  derivatives  with  respect  to  x;  and  r  t  in  Eq.  (D.2). 
Thus,  in  the  x  direction  we  have 


ISf  =  hr-  [<=03h- ‘1^3  Fl(0)B(cosh->^-,  firl  ) 

0 

-  J  ~~  Fm(xi  -  ^r,cosh<7)B((71(8ri  )d(7  (D.5) 

cosh” 1  ^ — 


If  the  condition  that  the  source  strength  at  the  body  apex 
is  equal  to  zero  is  now  imposed,  Fo(0)=0,  Eq.  (D.5)  is  then 
expressed  as 


0 

§£*-  =  "f  §7“  "  /5ricosh<r)B(71^ri)da 

i  J  v  i 

cosh*  1  — — 

fir. 


ID.  6) 


Eq.  (D.6)  after  discretization  along  the  body  axis  and  upon 
applying  the  Harmonic  Gradient  model  can  be  expressed  as 
i 


a*i  -  Y  f 

37)  ‘  Z  *4  e  J 


j  ♦  1 


cosCX/ffrjSinhcrlda 


(D.7) 


j  =  l 


Since  the  integrands  of  Eq.  (D.7)  are  regular  the  integrals 
of  Eq.  (D.7)  are  evaluated  numerically  by  a  Gaussian  gradrature. 
The  velocity  in  the  r  direction  is  given  by 

!fj-  =  a?7Ccosh'1^7]  Fa  (O)B(cosh-i^-,  Ar,  ) 


0 

-[  |^— [  Fm  (x  j  -/Sr  j  coshv)B  ( <r1  fir  i  )  ]  da  (D.8) 

cosh”1^- 

fir, 

Now  applying  again  Fn(0)  =  0  and  after  discretization  and 
making  use  of  the  Harmonic  Gradient  model  Eq.  (D.8)  can  be 
expressed  as 
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EXPLICIT  EXPRESSIONS  OF  MATRICES  [A],,  [Y]..  AND  [H]t 

IN  EQ.  (7.11) 
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The  explicit  expressions  of  matrices  [ A ] t ,  [Y]jl  and  [H]^,  involving 
all  the  elements,  in  Eq.  (7.11)  are  as  follows: 
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Fig.  2 


Sequence  of  flutter  modes  at  intervals  of  1/100 
second  during  the  flutter  incident. 
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TOTAL  VELOCITY/ FREESTREAM  VELOCITY 


Fig.  5  Mean  flow  pressures  for  a  parabolic-ogive  a 
M^-2.0  and  angle  of  attack  a=0* . 
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Mean  flow  pressures  for  a  parabolic-ogive 
M^=3.0  and  angle  of  attack  a  =  0*  . 
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Harmonic  gradient  model. 
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d  out-of-phase  pressure  coefficients 
in  pitching  mode  at  M  -2.0  und  reduced 


-phase  un<l  out-of-phase  pressure  coefficients 
r  u  cone  in  first  bending  mode  at  M^-2.0  and 

duced  frequency  k=2.0. 


1  36 


In-phuse  und  out-of-phase  pressure  coefficients 
for  a  cone  in  second  bending  mode  at  M  =2.0  and 


Modulus  and  argument  of  the  generalized  forces 
versus  Mach  number  for  a  cone  in  first  bending 
mode  (1=3)  at  reduced  frequency  k=2.0. 


PRESENT  HPP 
LINEARIZED  MOC 
SLENDER  BODY 


Comparison  of  theoretical  damping-in-pitch  moment 
coefficients  for  a  parabolic  ogive  at  various  Mach 
□umbers . 
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MACH  NUMBER,  M  <*) 


Fig.  20  Comparison  of  theoretical  and  experimental 
damping- in-pi tch  moment  coefficients  for  a 
parabolic-ogive  cylinder  at  various  Mach  numbers. 
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Comparison  of  theoretical  and  experimental 
damping-in-pi tch  moment  coefficients  for  a  cone 
frustrum  at  various  Mach  numbers. 


Fig.  22  Modulus  uud  urgument  of  the  generalized  forces 
versus  reduced  frequency  for  a  cone  and  a 
parabolic  ogive  in  pitching  mode  at  M  =2.5. 
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Pressure  coefficients  for  a  Saturn  SA-1 
configuration  at  M^=2.0,  reduced  frequency  k=1.8 

and  center  of  rotation  at  the  apex. 


53 


Comparison  of  computed  and  experimental 
ucrodynumic  dumping  coefficients  versus  Much 
number  for  a  Saturn  SA-1  configuration  vibruting 
in  first  bending  mode. 


u CENTER  OF 
GRAVITY  OF 
MODEL 


FLUTTER  SPEED,  V, /( b/2  cuayu) 


Flutter  speed  boundaries  versus  Mach  number  for  a 
7.5*  cone  . 
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FLUTTER  FREQUENCY, 


Fig.  34 


Flutter  frequency  boundaries  versus  Mach  number 
for  a  7.5*  cone . 


MACH  NUMBER,  M^ 


Fig.  35  Divergence  boundaries  versus  Mach  number  for  a 
7.5*  cone  . 
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8  (DEG.) 


Fig.  :57 


Pressure  coefficient  for  an  asymmetric  cone  at 
M  =2.0,  a  =  0  deg. 
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8  (DEG.) 


Pressure  coefficient  for  a  half  circle 
ellipse  cone  at  M  =2.0,  a=Q  deg. 


a/ b 


39 


Normal-force  coefficient  for  elliptic  cones  versus 
the  ellipticity  ratio  a/b  at  M  =3.0. 


D amp i n g- in-pi t ch  normal  force  coefficient  and 
damping-in-pitch  moment  coefficient  for  elliptic 
cones  versus  the  ellipticity  ratio  a/b  at  M  =3.0 
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Fig.  4  3  In-phase  and  aut-of-pliose  pressure  coefficient  for 
elliptic  cones  in  pitching  mode  ut  M  =3.0 

M 

and  reduced  frequency  k=1.0  along  the  xuxis. 


In  phase  and  out-of-phase  pressure  coefficient  for 
elliptic  cones  in  second  bending  mode  at  M  =3.0 

O* 

and  reduced  frequency  k=1.0  along  the  x-axis  • 
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liuse  and  out-of-phuse  pressure  coefficient  for 
E>  t  i  c  cones  in  second  bending  mode  at  M  =3.0 

•m 

reduced  frequency  k=l.O  along  the  polar 


Fig.  48  Pressure  distribution  per  unit  angle  of  attack 
along  the  y-axis  for  a  circular  cone,  elliptic 
cone  and  a  flat  plate  at  M  =2.0. 
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BUNDLED  TRIPLET  (BTM)  METHOD 
FOR  SUPERSONIC  FLOW 


owing  uunj1 ed  Triplet  Arrangements. 


omain  of  Influence  Based  on  the 
od  . 


PRESENT  METHOD 
DEVAN  {REF.  35)  X=l 
USSAERO  (REF  2  ) 


Steady  Pressure  Distributions  for  an  Asymmetric  Cone 
M  =  2.0  and  Angle  of  Attack  a  =  0°. 
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Fig.  53  Steady  Pressure  Distributions  for  an  Elliptic  Cone  a 
M  =  2.0  and  Angle  of  Attack  a  =  5°  . 
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Real  Part  of  the  Generalized  Aerodynamic  Force  Qj  l 
Versus  Circumferential  Mode  Number. 


< 


a 


n 


Fig.  5  i  . 


Imaginary  Part  of  the  Generalized  Aerodynamic  Force  Q 
Versus  Circumferential  Mode  Number. 


:8 


Sketches  of  Wing-Body  Configurations 

a)  Aspect  Ratio  AR  =  3  Triangular  Wing-Body 

b)  Aspect  Ratio  AR  =  3  Swept  Wing-Body 
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Static  and  Dynamic  Moment  Derivatives  for  an  Aspect 
Ratio  AR  =  3  Triangular  Wine-Body. 


tatic  and  Dynamic  Moment  Derivatives  for  an  Aspect 


CONICAL  COORDINATES  s,  v,  9 
CYLINDRICAL  COORDINATES  X,  r,  9 


Fig.  61  conical  Coordinate  and  Cylindrical  Coordinate  Systems 
for  a  Circular  Cone. 
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